{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is the reproduction of an exercise found at http://people.ku.edu/~gbohling/cpe940/Kriging.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "sys.path.append('../geostatsmodels')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geostatsmodels import utilities, variograms, model, kriging, geoplot\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll read the data from `ZoneA.dat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = utilities.readGeoEAS('../data/ZoneA.dat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want the first, second and fourth columns of the data set, representing the x and y spatial coordinates, and the porosity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = z[:,[0,1,3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be interested in determining the porosity at a point (2000,4700)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pt = [2000, 4700]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot our region of interest as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VdW5//HPNwnIpAIySAiIWigoYEpTBbXU4YLKVVEvVrnYoth69dYrpSrqtVa0pba2lmKt+nNA9GqlFqtwvRQHJq2iGDQiiigKyqTMgwKBkOf3x16Jh5DhhJxzcnLyvF+v/crea09rZXjOytprryUzwznnXGbJqu8MOOecSzwP7s45l4E8uDvnXAby4O6ccxnIg7tzzmUgD+7OOZeBPLi7BkXSpZL+Wd/5aGgkdZX0paTs+s6LSw0P7g2UpBHhj7XiYpJ+UU95Ghfuf0INx7WWNEnS55K2S/pQ0o2pymcVeTJJ36hm/3GS3pO0QdLPYtKbSHpDUpdqzu0Wrl/2M1qR6vKa2Wdm1srM9oY8zZX0o1TmwaWWB/cGysyeCH+s5QvwU+AL4MFU50eSgB8Cm8LX6kwAWgG9gEOBc4FlSc1g3d0BXAccB9ws6fCQ/jPgaTNbGcc1Woef03DgF5LOrE0GFPG/WRcX/0XJEJK+BfwRuNjM1oa0XEnTJW2StEzSj2OOHyfpKUmPhdrze5IKYvbnSnpa0npJyyVdU0MWvgt0Aq4BLpbUtJpjvwP8xcw2m1mpmX1gZlPDfctquTkxealYy5SkeyRtlfSBpNNjdlwq6ZNQpuWSRsTsGyVpiaTNkp6XdERIfzkc8k6oWV9USZ6PBGab2WrgI6BrOP/fiD6s4mZm84H3gN7h/idKejOU501JJ1Yo+3hJrwI7gKNq+LkeL6lQ0jZJX0j6Q8Xvq6TxRD+ve0J575H0Z0l3xeYz3GNMbcrm0oiZ+dLAF6A18DFwQ4X0l4F7gWZAPrAeOC3sGwfsAoYA2UQ109fDvixgIfALoClwFPAJcEY1eXgYeApoAmwE/q2aYx8iCm6XAd0r7OsGGJATkzYX+FFYvxQoAcaEe10EbAXaAi2BbcA3w7GdgGPD+lCi/w56ATnAz4HXYu5hwDeqyfPfgHOAPOBz4DDgWeB7cfx8yssECDiJKFCfHvK9GfhB2D88bB8WU/bPgGPD/iY1/FznAz8I662A/pV9X2O/p2H7eGANkBW224U8dqzv329fDmyp9wz4UscfYBQspgPTAMWkdwH2AgfHpN0BTA7r44CXYvYdA+wM6ycAn1W4z03AI1XkoUUIqueF7f8HTKsmz82B/yb6ANkTgu5ZYd8+QSiklQciouC+pkJZF4Tg2BLYQlSbbl7hnv8ALo/ZzgrB64iwXVNwPwKYAbwVAvC5wP8AXcP3fh5wYRXnlpVpSwjcS4Brwr4fAAsqHD8fuDSm7LfX4uf6MnAb0K6KPFQa3EPaEmBQWL8amFHfv9++HPjizTIN3w1EtbqRFv4qg1xgk5ltj0n7FOgcs/15zPoOoFloDjkCyJW0pWwhCsYdq8jD+US16Rlh+wngLEntKzvYzHaa2a/N7NtENeCngL9JahtHeQFWVyjrp0CumX1FVJO/Elgr6f8k9QzHHAFMjCnPJqIPxtjvR5XM7FMzG2Jm/YiC+S+J2uB/D/yVKNj/oYYytDOzNmbWy8zuDmm5If+xKv6cYtvza/q5Xg70AD4ITTxnx1O+4FHgkrB+CdGHl2ugPLg3YJJOAW4GhpnZlgq71wBtJR0ck9YVWB3HpVcCy82sdcxysJkNqeL4kURNAJ9J+pyoCaMJ8O813cjMtgG/Jqp1Hwl8FXa1iDns8AqndQ4PcMt0JSovZva8mQ0iapL5gK8fLq8E/qNCmZqb2Ws15bESvwAeNLMvgD5AoZltBVYBVfa4qcIaog+eWBV/TrEfZNX+XM3sIzMbDnQAfgtMldSykvtWNhzs48BQSccRNV89W5uCuPTiwb2BktQJmAL81Mzerrjfot4brwF3SGomqS9Rre7xOC6/ANgu6QZJzSVlS+ot6TuV5KMzUdvx2UTtv/lEPUp+SxW9ZiTdIuk7kppKagaMJmqyWGpm64kC1SXhvqOAoytcogNwTeiGeCFRIJohqaOkoSGYFQNfAqXhnPuBmyQdG/JwaDi3zBdEzxaqJekY4BTgvpC0HDhNUkegO1H7eG3MAHpI+vfwsPMioiay5yo7uKafq6RLJLU3s1Ki7yl8/T2ItV95zWwV8CZRjf1pM9tZy7K4dFLf7UK+HNhCVHs0ogBWcbk/HJNHFCQ2ET1wvTLm/HHA4zHb3di3TTYXeJKo6WYz8DrwL5Xk40ZgYSXpuUTt6b0r2fdzYDFRO/0movbfE2P2n0UUNLcAdxG1Z8e2ub8K3EP0IPVDYHDY1ykcuzWcOxc4Jua6PwDeDfddCUyK2XclsDac9/1qvu9zgBNito8D3gc2AD+r4px9vreV7D+Z6PnD1vD15Jh9c9m/bby6n+vjwLrwe/AeXz8HqfjzHRC+d5uBu2POvyQcd2p9/477UrdF4QfqnHNIGkj0AXGEeXBo0LxZxjkHRG/bEjWRPeSBveHz4O6cQ1IvoiapTkQvw7kGzptlnHMuA3nN3TnnMlBOzYc0PO3atbNu3brVdzaccw3AwoULN5hZpS/cxUtSbZpAnjezWg0adyAyMrh369aNwsLC+s6Gc64BkFTxDeFka5eKm2RkcHfOuVRTnOHUKElyTiIe3J1zrs6yEM3iOtL4Msl5iXhwd865OhPpFk7TKzfOOdcgKe5mmVRJr9w451yDlV7hNL1y45xzDZLX3J1zLuPIg7tzzmWi+HvLpIoHd+ecS4j0CqfplRvnnGuQvFnGOecykPdzd865DOQ1d+ecy0BKuweqSR3PXVJrSVMlfSBpiaQBktpKelHSR+Frm3CsJN0taZmkRZL6xVxnZDj+I0kjk5ln55yrrbKukPEsqZLsyTomAjPNrCfRLPFLgBuBWWbWHZgVtiGa8b57WK4A7gOQ1Ba4FTgBOB64tewDwTnn0kNZm3s8S2okLbhLOhQYCDwMYGa7zWwLMBR4NBz2KHBeWB8KPGaR14HWkjoBZwAvmtkmM9sMvAgkfaB755yLX+OquR8JrAcekfS2pIcktQQ6mtnacMznQMew3hlYGXP+qpBWVfo+JF0hqVBS4fr16xNcFOecq04jqrkTlaIfcJ+ZfQv4iq+bYACwaHbuhMzQbWYPmFmBmRW0b1+nGbOcc66WGlfNfRWwyszeCNtTiYL9F6G5hfB1Xdi/GugSc35eSKsq3Tnn0kTUWyaeJVWSFtzN7HNgpaRvhqTTgfeB6UBZj5eRwLSwPh34Yeg10x/YGppvngcGS2oTHqQODmnOOZcm0q/mnuw7/RfwhKSmwCfAZUQfKE9Juhz4FPh+OHYGMARYBuwIx2JmmyT9EngzHHe7mW1Kcr6dc64WGtkbqmZWBBRUsuv0So414CdVXGcSMCmxuXPOucTwIX+dcy5jpVc4Ta/cOOdcg+TjuTvnXAbyZhnnnMtAjeyBqnPONRZec3fOuYzjzTLOOZeBvFnGOecyjtKwt0yyx3N3LiONGjWKDh060Lt37/K066+/np49e9K3b1/OP/98tmzZUum5W7ZsYdiwYfTs2ZNevXoxf/58AIqKiujfvz/5+fkUFBSwYMGClJTFJUa6DT/gwd25A3DppZcyc+bMfdIGDRrE4sWLWbRoET169OCOO+6o9NzRo0dz5pln8sEHH/DOO+/Qq1cvAMaOHcutt95KUVERt99+O2PHjk16OVyiNK4hf53LWAMHDqRt27b7pA0ePJicnOiPt3///qxatWq/87Zu3crLL7/M5ZdfDkDTpk1p3bo1AJLYtm1b+XG5ubnJLIJLqMY3cJhzjdKkSZO46KKL9ktfvnw57du357LLLuOdd97h29/+NhMnTqRly5b88Y9/5IwzzuC6666jtLSU1157rR5y7g5M+vWW8Zq7czFKthSx+6PfU7zkFnZ/ci+lO1bU+hrjx48nJyeHESNG7H/9khLeeustrrrqKt5++21atmzJb37zGwDuu+8+JkyYwMqVK5kwYUJ57d41BN4s41zaKtk0n71rpmK710Ppbmznp+xZ8RClOz6N+xqTJ0/mueee44knnkDSfvvz8vLIy8vjhBNOAGDYsGG89dZbADz66KNccMEFAFx44YX+QLVBaUSTdTjXkJiVsveL58H2VNixh5IvZlZ+UgUzZ87kzjvvZPr06bRo0aLSYw4//HC6dOnC0qVLAZg1axbHHHMMALm5ucybNw+A2bNn07179wMsjUs1eZu7c2lq71f7B/bAitfulzZ8+HDmzp3Lhg0byMvL47bbbuOOO+6guLiYQYMGAdFD1fvvv581a9bwox/9iBkzZgDwpz/9iREjRrB7926OOuooHnnkEQAefPBBRo8eTUlJCc2aNeOBBx5IUmFd4qXfS0yK5sjILAUFBVZYWFjf2XANiJWWsPuD28B277dPzXJpevToesiVSwVJC82sskmF4tZUPaw9f4rr2DWcWef7xcObZZwDlJVD9mEngZpU2NGE7A6D6ydTrgFJXLOMpEmS1klaXCH9vyR9IOk9SXfWdJ30+j/CuXqU3WEwKIu9G/8JpXsgpyXZHYeQfXCv+s6aS3tK5MPSycA9wGPlV5dOBYYCx5lZsaQONV3Eg7tzgZRFTofBZLf/FyjdDVkHVdrjxbn9Ja7N3cxeltStQvJVwG/MrDgcs66m63izjHMVSFkou5kHdlcLtWqWaSepMGa5Io4b9AC+K+kNSfMkfaemE7zm7pxzdVarmvuGA3igmgO0BfoD3wGeknSUVdMjxoO7c87VkUj6TEyrgL+HYL5AUinQDlhf1QneLOOcc3WW9OEHngVOBZDUA2gKbKjuBK+5O+dcnWWRlaDeMpKeBE4haptfBdwKTAImhe6Ru4GR1TXJgAd355xLgIT2lhlexa5LanMdD+7OOZcAjWrIX0krJL0rqUhSYUgbJ2l1SCuSNCTm+JskLZO0VNIZMelnhrRlkm5MZp6dc6720m/I31Tc6VQzq9jwP8HMfh+bIOkY4GLgWCAXeCk8OAD4MzCI6Inxm5Kmm9n7Sc63c87FKf0m60in3AwFpoQ3sJZLWgYcH/YtM7NPACRNCcd6cHfOpYn0GxUy2V0hDXhB0sIKb2FdLWlRGCCnTUjrDKyMOWZVSKsqfR+Srih742v9+iq7fjrnXMIp9JaJZ0mVZAf3k82sH3AW8BNJA4H7gKOBfGAtcFcibmRmD5hZgZkVtG/fPhGXdM65uInsuJZUSer/EWa2OnxdJ+kZ4Hgze7lsv6QHgefC5mqgS8zpeSGNatKdc67eCchWnHNjpGgKjaTV3CW1lHRw2TowGFgsqVPMYecDZWMWTwculnSQpCOB7sAC4E2gu6QjJTUleug6PVn5ds65A5Eti2tJlWTW3DsCz4SR9XKAv5jZTEn/Iymf6PNrBfAfAGb2nqSniB6UlgA/MbO9AJKuBp4HsoFJZvZeEvPtnHO1Uquae4okLbiH3i3HVZL+g2rOGQ+MryR9BjAjoRl0zrkEkYymWaXxHbw3uXkpk159d5xzroHKaiw1d+ecaywEKewHEx8P7s45V0eNqs3dOecaEw/uzjmXYURquznGw4O7c87VkUT8vWVSxIO7c87Vkbe5O+dchvLg7pxzGcbb3J1zLgNJXnN3zrmM5C8xOedchhHeW8Y55zKOt7k751wG8q6QzjmXidLwgWqy51B19WTUqFF06NCB3r1777fvrrvuQhIbNmyo8vxt27aRl5fH1VdfXZ52880306VLF1q1apWUPDtXprLf3+uvv56ePXvSt29fzj//fLZs2VLpuRMmTODYY4+ld+/eDB8+nF27dgHw3e9+l/z8fPLz88nNzeW8885LWH7Lau7pNBOTB/cMdemllzJz5sz90leuXMkLL7xA165dqz3/lltuYeDAgfuknXPOOSxYsCCh+XSuMpX9/g4aNIjFixezaNEievTowR133LHfeatXr+buu++msLCQxYsXs3fvXqZMmQLAK6+8QlFREUVFRQwYMIALLrggYfkV0WQd8Syp4sE9Qw0cOJC2bdvulz5mzBjuvPNOwvSHlVq4cCFffPEFgwcP3ie9f//+dOrUqYqznEucyn5/Bw8eTE5O1JLcv39/Vq1aVem5JSUl7Ny5k5KSEnbs2EFubu4++7dt28bs2bMTX3OPc0kVD+4N0N4tuylZuxMrrd2/eNOmTaNz584cd9x+sx+WKy0t5dprr+X3v/99XbPpXKWstITSnRuxkl0HfI1JkyZx1lln7ZfeuXNnrrvuOrp27UqnTp049NBD96ukPPvss5x++ukccsghB3z/yqRbs0xcD1QltQFygZ3ACjNLrw6djcTezbvZ+Ov32f3+NpQl1CqHttd9k2bf3r+GXtGOHTv49a9/zQsvvFDtcffeey9DhgwhLy8vUdl2rtyeNa9T8tksMAMzstv3psnR56Cs+Pt2jB8/npycHEaMGLHfvs2bNzNt2jSWL19O69atufDCC3n88ce55JJLyo958skn+dGPfpSQ8pRpUG+oSjoU+AkwHGgKrAeaAR0lvQ7ca2ZzUpJLh5mx/oZFlKz8CvaCYVjxbjaOe48O932bJnktqj3/448/Zvny5eW19lWrVtGvXz8WLFjA4YcfXn7c/PnzeeWVV7j33nv58ssv2b17N61ateI3v/lNUsvnMt/eje9T8ulLULrn67QNi0FZNP3G0LiuMXnyZJ577jlmzZpVadPiSy+9xJFHHkn79u0BuOCCC3jttdfKg/uGDRtYsGABzzzzTAJKtK8GE9yBqcBjwHfNbJ/H0pK+DfxA0lFm9nAyM+giez7czt7Pd+43c7qVlPLl9NW0+c/u1Z7fp08f1q1bV77drVs3CgsLadeu3T7HPfHEE+XrkydPprCw0AO7S4g9K1/eJ7ADUFrC3nWLsCPPQtlNqz1/5syZ3HnnncybN48WLSqvzHTt2pXXX3+dHTt20Lx5c2bNmkVBQUH5/qlTp3L22WfTrFmzOpcnVjr2c6+yzd3MBpnZ/1QM7GHfQjP7qQf21Nm7YTdkVfIQdC+UrNm/7XL48OEMGDCApUuXkpeXx8MPV/2jKiwsjOvf1LFjx5KXl8eOHTvIy8tj3LhxtSmCa+Rs97bKd0hYyc59kir7/b366qvZvn07gwYNIj8/nyuvvBKANWvWMGTIEABOOOEEhg0bRr9+/ejTpw+lpaVcccUV5dedMmUKw4cPT3jZ0rG3jMxq/rSR1BfoRkxN38z+nrxs1U1BQYEVFhbWdzYSqmTdLj6/bAHs2ffnpYOyOGRkNw4e1qWecuZcfIqXPEnppqVAhZiT05xmx49Fqp/+HZIWmllBzUdWLe+gvnb14c/FdexNnx1R5/vFo8anGJImAX2B94Cyjx0D0ja4Z6KcDs1oMehwds76AisOP4YckXVwDi3P8u6JLv01OeJ0ird8EppmQoDPakKTbmfUW2BPpHRrlonnEXV/Mzsm6TlxNWpzTXea9mjFl8+uxnbspdmJh3HIvx9BVksfRcKlv6wWHTjouCvY89kcbPtKOOhQmnT5Htltqn9e1BAIyGqAwX2+pGPM7P3aXlzSCmA70WPAEjMrkNQW+CtRM88K4PtmtlnRo++JwBBgB3Cpmb0VrjMS+Hm47K/M7NHa5iUTKEu0GpJLqyG5NR/sXBrKatGeg3p+v76zkRSJekEptJacDawzs94hbRzwY6JeiwD/bWYzqrtOPP8LPUYU4JdKWiTpXUmLapHXU80sP6aN6UZglpl1B2aFbYCzgO5huQK4LxSqLXArcAJwPHBr6HfvnHNpISuxD1QnA2dWkj4hxNL8mgI7xFdzfxj4AfAuX7e518VQ4JSw/igwF7ghpD9m0RPe1yW1ltQpHPuimW0CkPQiUcGfTEBenHOu7hL4EpOZvSypW12vE0/Nfb2ZTTez5Wb2adkS5/UNeEHSQkll/ZE6mtnasP450DGsdwZWxpy7KqRVlb4PSVdIKpRUuH79+oq700ZdRmvMzs4uH9Xu3HPPLU+fPXs2/fr1o3fv3owcOZKSkpKk5d85t79ajgrZrixWheWKGi5f5urQejIpntaLeIL725L+Imm4pAvKljgzc7KZ9SNqcvmJpH2GGQy19IR83JnZA2ZWYGYFZW+npaO6jNbYvHnz8lHtpk+fDkRjwYwcOZIpU6awePFijjjiCB59tFE+knCuXtUiuG8oi1VheSCOy98HHA3kA2uBu2o6IZ7g3hwoBgYD54Tl7DjOw8xWh6/rgGeI2sy/CM0thK9lr02uBmI7a+eFtKrSG6S6jNZYmY0bN9K0aVN69OgBRMOiPv300wnJq3MuPmXT7CVr4DAz+8LM9oZxvR4kiqXVqjG4m9lllSyjajpPUktJB5etE304LAamAyPDYSOBaWF9OvBDRfoDW0PzzfPAYEltwr8ig0NaxohntEaAXbt2UVBQQP/+/Xn22WcBaNeuHSUlJZS9tDV16lRWrlxZ3WWccwmW7Mk6yirEwflEsbRa1Q0c9nOiwcE2VbH/NKCFmVX1WlZH4JlQE80B/mJmMyW9CTwl6XLgU6CsX9QMom6Qy4i6Ql4GYGabJP0SeDMcd3tVeapPe7eXsPnBVXw1ayNkQauz2tNmVGeymlXfQSre0RoBPv30Uzp37swnn3zCaaedRp8+fTj66KOZMmUKY8aMobi4mMGDB5OdncpRo51zEjRN0ANVSU8SdSRpJ2kVUW/BUyTlEzVjrwD+o6brVNdb5l3gfyXtAt7i61EhuxO1+7wE/Lqqk83sE2C/qqiZbQROryTdiEahrOxak4BJ1eS1XllJKWuueI89q3eVDw+w9a9r2fX2NnIfOLbappZ4R2uEaKxqgKOOOopTTjmFt99+m6OPPpoBAwbwyiuvAPDCCy/w4YcfJqOYzrkqJHLgMDOrbPCbWo/jVd3AYdPM7CTgSqKhB7KBbcDjwPFmNsbM0rdbSgp99fJmSr4o3nfcl93G7o93sOvt7dWeWzZa44oVK1ixYgV5eXm89dZb+wX2zZs3U1xcDETDlr766qscc0z04nDZaI/FxcX89re/LR9QyTmXOuk2WUc8be4fmdlkM7vDzP5oZs+b2c6azmtMit//Etu5/ysAtscoXvrVPmkHOlrjkiVLKCgo4LjjjuPUU0/lxhtvLA/uv/vd7+jVqxd9+/blnHPO4bTTTktg6ZxzNUn2A9UDylM8o0I2NKkeFXLbs1+wceKn2K59A7xaZNHhlm/Q8pSaZ0pyztWPRIwK2b35sXb3N/4a17FDFvdJyaiQDX8otjTQatBhqImihrcyWZDVMpsWJ7eut3w551KjbJq9dKq5e3BPgKyWOeQ+0Jum32wJOYIc0azPweT+v94ox7/FzmU6QdpN1hHPeO53V5K8FSg0s2mV7GuUmnZrTt4jfdi7rQRlQVYrH4bXucairM09ncQTgZoBPYG/he1/A5YDx0k61cx+mqzMNUTZh3hQd64xaojBvS9wkpntBZB0H/AKcDJRX3jnnGvUROLGc0+UeIJ7G6AVUVMMQEugrZntlVSctJw551wDoQQO+Zso8QT3O4EiSXOJPqAGAr8O48W8lMS8OZcwH3zwAaNGjWL79u20bduWp59+mnbt2tV3tlyGKHugmk7ieYnpYeBE4FmikR1PNrOHzOwrM7s+2Rl0LlEef/xx3n33XU488UTuv//++s6Oyyjp9xJTvE//sojGlskBviHpG2b2cvKy5Vxi9ezZs3y9uLiYww47rB5z4zJNIseWSZR4ukL+FriIaHyZsv87DPDg7tKCme0zOFvF7VjPP/88//jHP5g/f36qsucagYba5n4e8E0z84enLu2MGzeOLVu2MGHCBCRhZowZM4bWrVszbty4fY4tLS3l8ssvZ86cObRu7W8Ou8RKt+Aez+uTnwBNkp0R52rLzNiyZQsTJ05kzJgx5YF94sSJbNmyhYrjJq1Zs4ZDDz2U7t2711OOXaZK9mQdByKemvsOot4ys4im2wPAzK5JWq6ci4MkJkyYAMDEiROZOHEiAKNHjy6vycdq06YNd91V49STzh0AQ1kNr+Y+Hfgl8BqwMGZxrt7FBvgylQV2gK1bt/LQQw+lKmuuMREoy+JaUqXGmruZPZqKjDh3IMqaYmKNGTOm0gCfm5vL1KlTU5k915g0lDZ3SU+Fr+9KWlRxSV0WnatcbBv76NGjKS0tZfTo0fu0wTuXMlkW35Ii1dXcR4evZ6ciI87VliRat269Txt7WRNN69atq5271rmEEmk3gHqVwd3M1oavn6YuO87Vzrhx4/bp114W4D2wu1RTQ2mWKSPpAkkfSdoqaZuk7ZK2pSJzzsWjYiD3wO5STgbZpfEtKRLvwGHnmNmSZGfGOecarDTrChlPcP/CA7tzztUgzf5hrDK4S7ogrBZK+ivRqJCxLzH9Pcl5c865BkGhn3s6qa7mfk7M+g5gcMy2AR7cnXOuTJo9UK2ut8xlAJJOMrNXY/dJOinZGXPOuQZDBjkNbLIO4E9xplVKUraktyU9F7YnS1ouqSgs+SFdku6WtCy8KNUv5hojQ4+djySNjPfezjmXKlJ8S6pU1+Y+gGgGpvaSfhaz6xBqNxfsaGBJOK/M9WZW8T3ws4DuYTkBuA84QVJb4FaggKg5aKGk6Wa2uRZ5cM655EqzNvfqau5NiSbGzgEOjlm2AcPiubikPOBfgXhGaxoKPGaR14HWkjoBZwAvmtmmENBfBM6M5/7OOZcSImqaiWdJkera3OdJ+ifQ18xuO8Dr/xEYS/ShEGu8pF8As4Abw0QgnYGVMcesCmlVpe9D0hXAFQBdu3Y9wOw659wBakA1d8xsL5B7IBeWdDawzswqDg98E9AT+A7QFrjhQK5fkZk9YGYFZlbQvn37RFzSOefiZEjxLakSz0tMRZKmA38DvipLjKOf+0nAuZKGAM2AQyQ9bmaXhP3Fkh4Brgvbq4EuMefnhbTVwCkV0ufGkW/nnEsNAdkNqOYeNAM2AqcR9X0/hzhGijSzm8wsz8y6ARcDs83sktCOjqIBQM4DFodTpgM/DL3TBUPxAAATqklEQVRm+gNbw+BlzwODJbWR1Iaov/3ztSmkc84lXQMa8hf4ur97Aj0hqT3RZ10RcGVInwEMAZYRvTR1Wbj/Jkm/BN4Mx91uZpsSnCfnnKuTdBsVssbgHnq8/ImomQXgFWC0ma2K9yZmNpfQlGJmp1VxjAE/qWLfJGBSvPdzzrmUEg3rgWrwCFGTSW5Y/jekOeecK5MV51IDSZMkrZO0uJJ910oySe3iyU5N2pvZI2ZWEpbJgHdHcS7JRo0aRYcOHejdu3d52vXXX0/Pnj3p27cv559/Plu2bIn7XIgmN+ncuTP5+fnk5+czY8aMpJah8Yizj3t8TTeTqeRdHkldiJ45fhbPReIJ7hslXRKGEciWdAnRA1bnXBJdeumlzJw5c5+0QYMGsXjxYhYtWkSPHj2444474j63zJgxYygqKqKoqIghQ4YkPN+NkkDZpXEtNTGzl4HKnitOIHpvKK5PiHiC+yjg+8DnwFqit1MT/ZDVOVfBwIEDadu27T5pgwcPJicnelTWv39/Vq2q/NFXZee6JIu/t0w7SYUxyxU1XVrSUGC1mb0Tb3bi6S3zKXBuvBd0zsVjL7CVaISPpgd0hUmTJnHRRRfV+rx77rmHxx57jIKCAu666y7atGlzQPd3MURtJuvYYGYFcV9aagH8N/sOu16j6gYO+0U155mZ/bI2N3LOlZkDTAVKiP7D/h7RqyDxj8c3fvx4cnJyGDFiRK3ufNVVV3HLLbcgiVtuuYVrr72WSZO8I1pCJK+3zNHAkcA7YX7gPOAtSceb2edVnVRdzf2rStJaApcDhwEe3J2rtUJgCrA7Ju1lomrfv8d1hcmTJ/Pcc88xa9asWk8G3rFjx/L1H//4x5x9do3vI7o4iOQNLWBm7wIdyu8lrQAKzGxDdedV2eZuZneVLcADQHOitvYpwFGJyLRzjc909g3shO15wJ4az545cyZ33nkn06dPp0WLFrW++9q1a8vXn3nmmf1607gDJKLJOuJZarqU9CQwH/impFWSLj+QLFX7QFVSW0m/AhYR1fL7mdkNZrbuQG7mnKvq5WojejH7a8OHD2fAgAEsXbqUvLw8Hn74Ya6++mq2b9/OoEGDyM/P58oroxe816xZs0/Pl8rOBRg7dix9+vShb9++zJkzhwkTJiSjkI2T4lxqYGbDzayTmTUJQ7g8XGF/t5pq7QCKXgytZIf0O+AColr7n83sy5qzlR4KCgqssLCwvrPhXCV+D7xXSXorYCLxdWBziSRpYW0ecFYm//Cj7cURlXdLrajDHy6q8/3iUd1v0rVEb6T+HFgjaVtYtkvaluyMOZeZhrF/75imRL2NPbA3aA1osg7/TXMu4boR9Wp7GlhB1DdhKJBff1lyCZDaER/jEc947s65hDoC+FmNR7kGpqGNCumcc64GAqXZZB0e3J1zLhG8WcY55zJMGo7n7sHdOefqLLU9YeLhwd055xJAada/0IO7c87VlYA4xmpPJQ/uzjmXCN4s45xzGUaG/IGqc85loNqNvpx0Htydcy4RvObunHMZxvu5O+dc5hEg7y3jnHMZyNvcnXMuwyj9hvxN+jtVkrIlvS3pubB9pKQ3JC2T9FdJTUP6QWF7WdjfLeYaN4X0pZLOSHaenXOutiSLa0mVVLwwOxpYErP9W2CCmX0D2AyUTf56ObA5pE8IxyHpGOBi4FjgTOBeSdkpyLdzzsUvy+JbUpWdZF5cUh7wr8BDYVvAacDUcMijwHlhfWjYJuw/PRw/FJhiZsVmthxYBhyfzHw751ytiLSbZi/ZNfc/AmOBssfIhwFbzKwkbK8COof1zsBKgLB/azi+PL2Sc8pJukJSoaTC9evXJ7oczjlXvRyLb0mRpAV3SWcD68xsYbLuEcvMHjCzAjMraN++fSpu6ZxzkTjb21PZ5p7M3jInAedKGgI0Aw4BJgKtJeWE2nkesDocvxroAqySlAMcCmyMSS8Te45zzqWHxtJbxsxuMrM8M+tG9EB0tpmNAOYAw8JhI4FpYX162Cbsn21mFtIvDr1pjgS6AwuSlW/nnDsgadbmXh/93G8Apkj6FfA28HBIfxj4H0nLgE1EHwiY2XuSngLeB0qAn5jZ3tRn2znnqiBS0/ewFlIS3M1sLjA3rH9CJb1dzGwXcGEV548Hxicvh845Vzc+/IBzzmWaFDe5xMODu3POJUKaPVD14O6cc3UkQD5wmHPOZRgfz9055zKUB3fnnMswMu8t45xzGcnb3J1zLgOlWbNMmr1T5RJl1KhRdOjQgd69e5en3XLLLfTt25f8/HwGDx7MmjVrKj137NixHHvssfTq1YtrrrmGaBSIr5177rn7XNc5R9oNP+DBPUNdeumlzJw5c5+066+/nkWLFlFUVMTZZ5/N7bffvt95r732Gq+++iqLFi1i8eLFvPnmm8ybN698/9///ndatWqV9Pw716CU9ZZJwGQdkiZJWidpcUzaLyUtklQk6QVJuTVdx4N7hho4cCBt27bdJ+2QQw4pX//qq69QJR1zJbFr1y52795NcXExe/bsoWPHjgB8+eWX/OEPf+DnP/95cjPvXAOUwCF/JxPNOhfrd2bW18zygeeAX9R0EW9zb2Bsbyl7PlmH7dpDk6M6kNXyoFqdf/PNN/PYY49x6KGHMmfOnP32DxgwgFNPPZVOnTphZlx99dX06tULiJp1rr32Wlq0aJGQsjiXMZS4iTjM7OXYOaRD2raYzZZAjTfzmnsDUrJmMxt/NY1tj/2T7U+9wcbx09gx74NaXWP8+PGsXLmSESNGcM899+y3f9myZSxZsoRVq1axevVqZs+ezSuvvEJRUREff/wx559/fqKK41xmib/NvV3ZrHFhuSKuy0vjJa0ERhBHzd2DewNhe0vZ8tBc7KtirLgEKy6BklK+euFd9ny6odbXGzFiBE8//fR+6c888wz9+/enVatWtGrVirPOOov58+czf/58CgsL6datGyeffDIffvghp5xySgJK5lxmUJbFtQAbymaNC8sD8VzfzG42sy7AE8DVNR3vwb2B2LN8PeypZBj7kr3sfH1ZXNf46KOPytenTZtGz5499zuma9euzJs3j5KSEvbs2cO8efPo1asXV111FWvWrGHFihX885//pEePHsydO/dAi+NcZkntBNlPAP9W00He5t5AWHFJFTvAdu7eL3n48OHMnTuXDRs2kJeXx2233caMGTNYunQpWVlZHHHEEdx///0AFBYWcv/99/PQQw8xbNgwZs+eTZ8+fZDEmWeeyTnnnJPMojmXGZJYVZbU3czKamdDgRrbY1WxD3MmKCgosMLCwvrORkKV7tjNxl9Ng5IKtfemORx8QQHN+nWrl3w519BJWmhmBXW5Rr+eufbqA/8R17Etvjeu2vtJehI4BWgHfAHcCgwBvgmUAp8CV5pZtXNJe829gchq0ZSWZx/HV//3ThTgDWiaTU5uaw46rmt9Z885l6AXlMxseCXJD1eSVi0P7g1IixN70KRLO3a9sYzSHbs5qE8eB/XtirL90Ylz9S7Nhh/w4N7ANOnSliZd9puC1jlXnxrrBNnOOZfJopmYvObunHOZx5tlnHMuw8jAJ+twzrkM5JN1OOdc5pE3yzjnXIYpG34gjSSt846kZpIWSHpH0nuSbgvpkyUtD4POF0nKD+mSdLekZWFQ+n4x1xop6aOwjExWnp1z7oAlaLKORElmzb0YOM3MvpTUBPinpH+Efdeb2dQKx58FdA/LCcB9wAmS2hK9fltA9F7mQknTzWxzEvPunHO1kNrAHY+k1dwt8mXYbBKW6ko/FHgsnPc60FpSJ+AM4EUz2xQC+ovsP0uJc87VH4GyLa4lVZL6TpWkbElFwDqiAP1G2DU+NL1MkFQ2lVBnYGXM6atCWlXpFe91Rdng9+vXr094WZxzrlqNaYJsM9sb5vzLA46X1Bu4CegJfAdoC9yQoHs9UDb4ffv27RNxSeeci08CJ8hOlJSMhmBmW4A5wJlmtjY0vRQDjwBlA6WsBrrEnJYX0qpKd865tJHACbITIpm9ZdpLah3WmwODgA9COzqSBJwHLA6nTAd+GHrN9Ae2mtla4HlgsKQ2ktoAg0Oac86lj6w4lxRJZm+ZTsCjkrKJivSUmT0nabak9kT/yBQBV4bjZxANSL8M2AFcBmBmmyT9EngzHHe7mW1KYr6dc652GtPwA2a2CPhWJemnVXG8AT+pYt8kYFJCM+iccwnko0I651wmSrN+7h7cnXOursp6y6QRD+7OOZcIPiqkc85lmvQbfsCDu3PO1ZEEaiy9ZZxzrlHx3jLOOZdhREpfUIqHB3fnnEsA7+funHMZxx+oOudc5vF+7s45l5lSORFHPDy4O+dcXaXhBNke3J1zLhG8WcY55zJNaqfQi4cHd+ecqyuBvObunHMZKM0eqCqaIyOzSNoOLK3vfNSTdsCG+s5EPWnMZYfGXf66lP0IM2tfl5tLmhnyEI8NZnZmXe4Xj0wN7oVmVlDf+agPXvbGWXZo3OVvzGWvSpqNhuCccy4RPLg751wGytTg/kB9Z6Aeedkbr8Zc/sZc9kplZJu7c841dplac3fOuUbNg7tzzmWgjAvuks6UtFTSMkk31nd+EkHSJEnrJC2OSWsr6UVJH4WvbUK6JN0dyr9IUr+Yc0aG4z+SNLI+ylJbkrpImiPpfUnvSRod0jO+/JKaSVog6Z1Q9ttC+pGS3ghl/KukpiH9oLC9LOzvFnOtm0L6Ukln1E+JakdStqS3JT0XthtFuRPGzDJmAbKBj4GjgKbAO8Ax9Z2vBJRrINAPWByTdidwY1i/EfhtWB8C/INonLr+wBshvS3wSfjaJqy3qe+yxVH2TkC/sH4w8CFwTGMofyhDq7DeBHgjlOkp4OKQfj9wVVj/T+D+sH4x8Newfkz4WzgIODL8jWTXd/niKP/PgL8Az4XtRlHuRC2ZVnM/HlhmZp+Y2W5gCjC0nvNUZ2b2MrCpQvJQ4NGw/ihwXkz6YxZ5HWgtqRNwBvCimW0ys83Ai0DS35KrKzNba2ZvhfXtwBKgM42g/KEMX4bNJmEx4DRgakivWPay78lU4HRJCulTzKzYzJYDy4j+VtKWpDzgX4GHwrZoBOVOpEwL7p2BlTHbq0JaJupoZmvD+udAx7Be1fegwX9vwr/b3yKqwTaK8oemiSJgHdEH0sfAFjMrCYfElqO8jGH/VuAwGmbZ/wiMBUrD9mE0jnInTKYF90bJov9BM7pPq6RWwNPAT81sW+y+TC6/me01s3wgj6jW2bOes5R0ks4G1pnZwvrOS0OWacF9NdAlZjsvpGWiL0JzA+HrupBe1fegwX5vJDUhCuxPmNnfQ3KjKT+AmW0B5gADiJqaykZ0jS1HeRnD/kOBjTS8sp8EnCtpBVHT6mnARDK/3AmVacH9TaB7eKrelOjhyvR6zlOyTAfKenyMBKbFpP8w9BrpD2wNzRfPA4MltQk9SwaHtLQW2k4fBpaY2R9idmV8+SW1l9Q6rDcHBhE9c5gDDAuHVSx72fdkGDA7/FczHbg49Co5EugOLEhNKWrPzG4yszwz60b0NzzbzEaQ4eVOuPp+opvohai3xIdEbZM313d+ElSmJ4G1wB6idsPLidoUZwEfAS8BbcOxAv4cyv8uUBBznVFED5WWAZfVd7niLPvJRE0ui4CisAxpDOUH+gJvh7IvBn4R0o8iClLLgL8BB4X0ZmF7Wdh/VMy1bg7fk6XAWfVdtlp8D07h694yjabciVh8+AHnnMtAmdYs45xzDg/uzjmXkTy4O+dcBvLg7pxzGciDu3POZSAP7q5eSdorqShmqfVInpJOkXRizPaVkn6YwDxOlXRULY7vI2lyou7v3IHIqfkQ55Jqp0Wv19fFKcCXwGsAZnZ/XTNVRtKxRCMJfhLvOWb2rqQ8SV3N7LNE5cW52vCau0tLkn4h6U1JiyU9EN5URdI1isZ2XyRpShhM7EpgTKj5f1fSOEnXhePnSvptGBf9Q0nfDektJD0VrvVMGAe8oJKsjODrNyGR9KWk3ykaX/0lSceHe3wi6dyY8/6X6O1K5+qFB3dX35pXaJa5KKTfY2bfMbPeQHPg7JB+I/AtM+sLXGlmK4jG9p5gZvlm9kol98gxs+OBnwK3hrT/BDab2THALcC3q8jfSUDsAFYtiV5vPxbYDvyKaFiA84HbY44rBL4b5/fAuYTzZhlX36pqljlV0ligBdEEG+8R1YYXAU9IehZ4Ns57lA02thDoFtZPJhqMCjNbLGlRFed2AtbHbO8GZob1d4FiM9sj6d2Ya0M0kFlunPlzLuG85u7SjqRmwL3AMDPrAzxINH4IRBM4/JloZqo3Y0YJrE5x+LqX2ldodsbcG2CPfT1mR2nZtc2stMK1m4VznasXHtxdOioLphvCOO7DACRlAV3MbA5wA9HQrq2ImkcOruU9XgW+H657DNCniuOWAN+o5bUBehAN9uVcvfBmGVffmoeZhsrMNLMbJT1IFBw/JxrKGaI5ch+XdCjR6I93m9kWSf8LTJU0FPivOO97L/CopPeBD4iafbZWctz/EfXGeamW5To1nOtcvfBRIV2jJCkbaGJmuyQdTRS8v2nR3LuxxzUnGkf8JDPbG+e1DwLmASfb19PCOZdSXnN3jVULYE6Y5UnAf1YM7ABmtlPSrURzb8bbZ70rcKMHdlefvObunHMZyB+oOudcBvLg7pxzGciDu3POZSAP7s45l4E8uDvnXAb6/2nl1RxUmXaiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(P[:,0], P[:,1], c=P[:,2], cmap=geoplot.YPcmap)\n",
    "plt.title('Zone A Subset % Porosity')\n",
    "plt.colorbar()\n",
    "xmin, xmax = 0, 4250\n",
    "ymin, ymax = 3200, 6250\n",
    "plt.xlim(xmin,xmax)\n",
    "plt.ylim(ymin,ymax)\n",
    "for i in range(len(P[:,2])):\n",
    "    x, y, por = P[i]\n",
    "    if (x < xmax) & (y > ymin) & (y < ymax):\n",
    "        plt.text( x+100, y, '{:4.2f}'.format( por ) ) \n",
    "plt.scatter(pt[0], pt[1], marker='x', c='k')\n",
    "plt.text(pt[0] + 100 , pt[1], '?')\n",
    "plt.xlabel('Easting (m)')\n",
    "plt.ylabel('Northing (m)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can determine the parameters for our model by looking at the semivariogram and trying to determine the appropriate range and sill."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 250\n",
    "lags = np.arange(tolerance, 10000, tolerance*2)\n",
    "sill = np.var(P[:,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The semivariogram plotting function, `svplot()`, plots sill as a dashed line, and the empirical semivariogram as determined from the data. It optionally plots a semivariance model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VNX5wPHvmz0BSVgtWyYoi4ILAooIEkBcoCziCqYqan9xrVLUKqVVq+LWutSKCtUqhQgBcQFEKQEUFxChoICIIBCCiIQdDBCSvL8/5iad7JNkJjPJvJ/nmYeZc8+9970zw7y555x7rqgqxhhjTHnCAh2AMcaY4GaJwhhjTIUsURhjjKmQJQpjjDEVskRhjDGmQpYojDHGVMgShQlpInKhiGysa9s2pjaJXUdhgpWI9AGeAboA+cAGYIyqfhXQwIwJMRGBDsCYsohII2AecDswE4gCLgSOBzIub4lIhKrm1bVtG1MWa3oywaojgKpOV9V8VT2qqv9R1W8ARORmEdkgIvtFZIGIuApXFBEVkTtEZJOIHBaRx0TkVBH5QkQOichMEYly6vYTkR3O8wdE5G3PIETk7yLyovP8Jmefh0Vki4jc6lGvn4jscLaxC3jDc9tOndNF5GMROSAi60VkmMeypiIy14nvKxF5XEQ+K3FMd4rIJmCTR2xZzjqrRORCj/qPiMgsEZnmxLtWRDqKyDgR2e2sd4lvPipT31miMMHqeyBfRKaIyCARaVy4QESGA38ErgCaA58C00usfynQHTgf+AMwGfgN0BY4AxhVxj5nAINF5CRnP+HANcBbzvLdwBCgEXAT8LyIdPNY/1dAE8AFpHpuWEQigbnAf4AWwO+ANBHp5FSZCPzibONG51HS5UBPoLPz+iugq7PPt4BZIhLjUX8oMBVoDKwGFuD+P98aeBSYVMY+jClNVe1hj6B8AKcDbwI7gDxgDnAy8CFwi0e9MCAHcDmvFejtsXwV8IDH62eBF5zn/YAdHss+A25wnl8M/FBBfO8B93hsJxeI8VhetG3czWa7gDCP5dOBR4Bw4ATQyWPZ48BnHq8VGFDJ+7UfONt5/giw0GPZUOAIEO68PsnZZkKgP2d7BP/DzihM0FLVDao6WlXb4D4LaAW8gPsv9r87TTgHgH2A4P5LudDPHs+PlvG6YTm7fYv/nW1cx//OJnDObJaLyD5nv4OBZh7rZqvqsXK22wrIUtUCj7JMJ+bmuPsLszyWeT4vs0xE7nOawg468cSXiKfkMe9R1XyP11D++2BMEUsUpk5Q1e9wn12cgfsH81ZVTfB4xKrqFz7Y1Sygn4i0AUbgJAoRiQZmA38DTlbVBGA+7gRVFGYF290JtBURz/9zicCPQDbuM6Y2HsvalrGNou07/RF/wN001tiJ52CJeIzxCUsUJiiJyGkicq/zg42ItMX9l/5y4FVgnIh0cZbFi8jVvtivqmYDHwNvAFtVdYOzKAqIxvlRF5FBQFU6g7/E3Tz2BxGJFJF+uJuDZjh/5b8DPCIicSJyGnBDJds7CXdyyQYiROQh3H0nxvicJQoTrA7j7rj9UkR+wZ0g1gH3quq7wNPADBE55JQP8uG+3wIG4tHspKqHgbtxD9Xdj7tZao63G1TVXNyJYRCwB3gZd1/Id06Vu3A3He3C3QE9nYqHAi8APsLd6Z8JHKPs5ipjaswuuDMmCInI08CvVLWs0U/G1Co7ozAmCDhNbWeJ23nALcC7gY7LGLArs40JFifhbm5qhXu00rPA+wGNyBiHNT0ZY4ypkDU9GWOMqVCda3pq1qyZJiUlBToMY4ypU1atWrVHVZtXZ906lyiSkpJYuXJloMMwxpg6RUQyq7uu35qeRORfziyV68pZLiLyoohsFpFvSkyuZowxJkj4s4/iTeCyCpYPAjo4j1TgFT/GYowxppr8lihUdSnuydrKMxz4t7otBxJEpKW/4jHGGFM9gRz11JriUw7soPjsn0VEJFVEVorIyuzs7FoJzhhjjFudGB6rqpNVtYeq9mjevFqd9sYYY6opkIniR4pPpdzGKTOmXklLSyMpKYmwsDCSkpJIS0sLdEjGVEkgE8Uc4AZn9NP5wEFV/SmA8Rjjc2lpaaSmppKZmYmqkpmZSWpqqiULU6f4bQoPEZmO+1aQzXDPXfMwEAmgqq+KiAAv4R4ZlQPcpKqVXiDRo0cPtesoTF2RlJREZmbp4esul4tt27bVfkAmZInIKlXtUZ11/XbBnaqWdfN6z+UK3Omv/RsTDLZv316lcmOCUZ3ozDamrkpMTKxSuTHByBKFMX40YcIE4uLiipXFxcUxYcKEAEVkTNVZojDGj1JSUrjzzuItrM8//zwpKSkBisiYqrNEYYyfHTlyhLi4OD788EMAWrVqFeCIjKkaSxTG+JGqMnfuXC655BKSk5OJiorik08+CXRYxlSJJQpj/GjNmjXs2LGDoUOHEhsby3nnncfSpUsDHZYxVWKJwhg/mjt3LiLCr3/9awCSk5NZtWoVhw8fDnBkxnjPEoUxfjR37lzOO+88Tj75ZAD69u1Lfn4+y5YtC3BkxnjPEoUxfrJz505WrlzJsGHDisouuOACwsPDrZ/C1CmWKIzxkw8++ACAoUOHFpU1bNiQ7t27Wz+FqVMsURjjJ3PmzMHlcnHGGWcUK09OTmbFihUcPXo0QJEZUzWWKIzxg5ycHDIyMhg6dCju+S//Jzk5mdzcXJYvXx6g6IypGksUxvjBokWLOHbsWLH+iUK9e/dGRKz5ydQZliiM8YO5c+dy0kknkZycXGpZQkICXbt2tQ5tU2dYojDGxwoKCpg3bx6XXnopUVFRZdbp27cvy5YtIzc3t5ajM6bqLFEY42P//e9/+emnn4qNdiopOTmZY8eO8dVXX9ViZMZUjyUKY3xszpw5hIWFMXjw4HLrXHjhhQDWT2HqBEsUxvjY3LlzueCCC2jWrFm5dZo1a0aXLl2sn8LUCZYojPGhrKws1qxZU2GzU6G+ffvy+eefk5eXVwuRGVN9liiM8aF58+YBeJUokpOTOXLkCKtXr/Z3WMbUiCUKY3xo7ty5nHrqqZx22mmV1u3bty+ANT+ZoGeJwhgfOXLkCIsWLSrzauyytGzZkg4dOliHtgl6liiM8ZGFCxeSm5tb5tXY5enbty+ffvop+fn5fozMmJqxRGGMj8ydO5f4+Hj69Onj9TrJyckcOHCAdevW+TEyY2rGEoUxPlBQUMAHH3zAoEGDiIyM9Hq9wik+rJ/CBDNLFMb4wIoVK9i9e7dXo508JSYm4nK5rJ/CBDVLFMb4wJw5cwgPD2fQoEFVXjc5OZmlS5eiqn6IzJias0RhjA/MnTuXCy+8kMaNG1d53b59+5Kdnc13333nh8iMqTlLFMbU0LZt21i3bl2Vm50KWT+FCXaWKKrho48+olOnTrRv356nnnqq1PLf//73dO3ala5du9KxY0cSEhIAWLJkSVF5165diYmJ4b333gPcN7rp1q0bXbt2pU+fPmzevLlWj8lU39y5cwHvrsYuy6mnnkrLli2tn8IEL1WtU4/u3btrIOXl5ekpp5yiP/zwgx4/flzPOussXb9+fbn1X3zxRb3ppptKle/du1cbN26sv/zyi6qqdujQQb/99ltVVZ04caLeeOONfonf+N7FF1+snTp1qtE2Ro4cqa1atdKCggIfRWVMccBKrebvrp1RVNGKFSto3749p5xyClFRUYwcOZL333+/3PrTp09n1KhRpcrffvttBg0aRFxcHAAiwqFDhwA4ePAgrVq18s8BhKC0tDSSkpIICwsjKSmJtLQ0n2370KFDfPzxx9U+myjUt29fdu7cyQ8//OCjyIzxnQh/blxELgP+DoQDr6nqUyWWJwJTgASnzoOqOt+fMdXUjz/+SNu2bYtet2nThi+//LLMupmZmWzdupUBAwaUWjZjxgzGjh1b9Pq1115j8ODBxMbG0qhRI5YvX+774ENQWloaqamp5OTkAO7PJDU1FYCUlJQab3/BggWcOHGiSldjl6Wwn2Lp0qW0b9++xnEZ40t+O6MQkXBgIjAI6AyMEpHOJar9CZipqucAI4GX/RVPIMyYMYOrrrqK8PDwYuU//fQTa9eu5dJLLy0qe/7555k/fz47duzgpptuKpZETPWNHz++KEkUysnJYfz48T7Z/ty5c2nSpAm9evWq0XZOP/10mjVrZh3aJij5s+npPGCzqm5R1VxgBjC8RB0FGjnP44GdfozHJ1q3bk1WVlbR6x07dtC6desy686YMaPMZqeZM2cyYsSIoit4s7Oz+frrr+nZsycA1157LV988YUfog8927dvr1J5VeTn5zN//nwGDx5MRETNTs5FhL59+/qtQ9ufzW91Yf+mhqrbuVHZA7gKd3NT4evrgZdK1GkJrAV2APuB7uVsKxVYCaxMTEz0cRdP1Zw4cULbtWunW7ZsKerMXrduXal6GzZsUJfLVWbnZM+ePXXx4sXFttm0aVPduHGjqqq+9tpresUVV/jvIEKIy+VS3H+QFHu0atWqxtv+9NNPFdD09HQfRKr6wgsvKKCZmZk+2V6hadOmaVxcXLHjj4uL02nTpvl0P8G6f+NGDTqzA50oxgL3Os97Ad8CYRVtN9CjnlRVP/jgA+3QoYOecsop+vjjj6uq6p///Gd9//33i+o8/PDD+sADD5Rad+vWrdqqVSvNz88vVv7OO+/oGWecoWeddZYmJyfrDz/84N+DqEXTpk1Tl8ulIqIul6tWfyCmTZumUVFRpRJF8+bNNTs7u0bbvv/++zUiIkIPHDjgk1hXr16tgE6dOtUn2ytUXrJ0uVw+3U+w7t+4BWui6AUs8Hg9DhhXos56oK3H6y1Ai4q2GwyJwngvGP6aPOecczQ8PLwoUT3yyCMaExOjvXv31mPHjlV7u6eddpoOHDjQZ3Hm5eVpQkKC/va3v/XZNlVVRaTMH2oR8el+ypKXl1fmvmtr/74SyD92fCVYE0WE88PfDogCvga6lKjzITDaeX467j4KqWi7lijqlkD/Nblv3z6NjIzU3//+98XK09PTFdCUlJRqXbuwadMmBfTvf/+7r0JVVdUhQ4Zox44dfbrNk08+uczPwN/NuIcPH9ahQ4eWmyjqyhlFMPyx4ws1SRR+Gx6rqnkichewAPfQ13+p6noRedQJeA5wL/BPEfm98wGMdg6oXBs3bqRfv37Fyq655hruuOMOcnJyGDx4cKl1Ro8ezejRo9mzZw9XXXVVqeW333471157LVlZWVx//fWllt97770MHTqUjRs3cuutt5Za/qc//YmBAweyZs0axowZU2r5E088wQUXXMAXX3zBH//4x1LLX3jhBbp27UpGRgaPP/54qeWTJk2iU6dOzJ07l2effbbU8qlTp9K2bVvS09N55ZVXSi1/++23adasGW+++SZvvvlmqeXz588nLi6Ol19+mZkzZ5Za/vHHHwPwt7/9reie0IViY2P58MMPAXjsscdYtGhRseWZmZmltldY3q9fP9q0acO0adMAGDNmDGvWrClWr2PHjkyePBmA1NRUvv/++2LLu3btygsvvADAb37zG3bs2FFseYMGDThx4gQjR47kyiuvZO/evUXL2rVrR1paGh07duShhx5i0KBBHD16tNj6Q4YM4b777gMo9r0r3M+BAwcAfPbdy8rKYsuWLVxwwQVERUXV+Lv3+OOPExsbW6ocYP/+/SxbtoxevXr5/LuXm5vL2rVr+eWXX7jxxhuZPn06ubm5RctFhNjYWHJycvz23WvatCmzZ88GYNy4cSxbtqzYcm+/e+WNnLv99tuLhliX9d3r1asXTz75JECp7x7ARRddxJ///GeAKn33CtX0d68q/HodhbqviZhfouwhj+ffAr39GYMJrLi4uFL/yQCio6NrZf8bNmygXbt2nHvuuaWWJSYmEh8fz8MPP1zlaxf27t1LXFwczZo181WoAEXTvRw4cIAWLVrUeHsffvgh27ZtY8SIEcyfP5/jx48THR1NfHw8u3fv5pZbbmHhwoU13o+nI0eOsHbtWvLy8khLS2PUqFFER0czZcoUjh8/TlhYGAUFBUUXmwYrVeX9998v94+dw4cP13JEAVTdU5FAPazpqW4pqzM5JiamVk7bs7OzNTw8vMxBBYWOHz+uycnJGhUVpZ999plX2923b5+Gh4fruHHjfBVqkRMnTmiDBg30jjvuqPG2Dh48qCeffLL26tWrzOa1OXPmaMOGDbVVq1a6atWqGu9PVXX+/PnasGFDbd26ta5Zs6bMOnv37tXWrVvraaedVjSFTTA5evSoTp48WTt16qSAhoeHB6TpztcIxj4Kfz0sUdQ93bp107CwsKJOVV92AFdk0qRJCuh///vfCuvt2bNHO3TooM2aNfNqtNlbb72lgH7xxRe+CrWYSy65RLt06VLj7TzwwAMK6IoVK8qt8/XXX2tiYqLGxsbq7Nmza7S/l19+WcPCwvScc87RHTt2VFh34cKFCuidd95Zo3360t69e/Xxxx8v6tPp1q2bTp8+Xf/973+X6qMAtGvXrnrw4MFAh+01SxQmaO3Zs0cjIyN17NixqqqakpKi8fHxmpOT4/d9DxgwQDt06OBVZ/X333+vTZo00dNOO033799fYd1Ro0Zp8+bNNS8vz1ehFjNhwgQFajR8d/PmzRoVFaU33HBDpXV37dql559/vgI6YcKEKnfu5+Xl6dixYxXQIUOG6OHDh71ab8yYMQro/Pnzq7Q/X9u6davefffd2qBBAwV00KBBumjRomLvg+eop8TERL3xxhs1PDxcTz/9dP3+++8DGL33LFGYoDVx4kQFdPXq1aqqumjRIgU0LS3Nr/vdtWuXhoWF6Z/+9Cev1/n44481MjJSL7roIs3NzS2zTm5uriYkJOjo0aN9FWophRfyvfPOO9XexogRI7RBgwb6448/elX/6NGjet111ymg119/vdfDho8cOaLDhw9XQO++++4qJc+jR49qly5d9Fe/+lWNr2mpTFnDW1euXKkjR47U8PBwjYyM1BtvvFHXrl3r9TYXL16sTZs21YSEBF2wYIEfo/fN8FxLFCZo9erVS88444yiv87y8/O1Xbt2OmDAAL/u96WXXlKgSv/xVVXffPNNBfS3v/1tmX9ZL1myRIEaN9NU5NixYxoTE6Njxoyp1vqFyXjChAlVWq+goEAfe+wxBfSCCy7Qn3/+ucL6O3fu1O7du2tYWFi1hwmvWbNGo6Ki9PLLL/fbFOtlDW8NCwtTQBs1aqT3339/pU1l5dmyZYueeeaZGhYWps8++6xfjsFXw3MtUZigVHitwdNPP12s/NFHH1VAt2zZ4rd99+nTp9rt/OPHj1dA//rXv5ZaNnbsWI2KivK6eaW6+vXrp+ecc06V1ztx4oSeeeaZmpSUpEePHq3WvmfOnKkxMTGalJRUbqL95ptvtG3bttqgQQOdO3dutfZT6JlnnlFAX3/99RptpzzlXcuTkJDgk6vqDx8+rCNGjFBAb7jhhmq/7+VJTEz0yXUolihMUHr44YdVRDQrK6tY+fbt21VE9KGHHvLLfrOyshTQRx99tFrr5+fn6zXXXKMiUqz5p6CgQNu3b6+XXXaZr0ItV+F7V9UfsldeeUUBnTVrVo32v2LFCm3ZsqWedNJJ+sEHHxRr+mjRooVGR0drq1atKh0o4I28vDzt16+fNmzYUDdv3lzj7ZVUG1em5+fn6yOPPKKA9uzZ0+smv4qsX79e7733Xp9d2W6JwgSdgoICPfXUU8ttYrr00ku1bdu2fukQfu655xQommSxOnJycrRnz54aGxurK1euVFX3RI+ATpw40Vehlmvx4sUK6Lx587xeZ9++fdq0aVNNTk72SRNIVlaWnnPOOQpoZGRkqR+pF198scb7KJSZmanx8fHaq1cvPXHihM+2u2zZsnKHt/rjyvDZs2drgwYNtGXLlrp8+fIqr3/gwAF99dVXtWfPngpoRESExsbG2hlFVR+WKOqGzz//XAF94403ylxeOIWGPzoBe/bsqV27dq3xdnbt2qUul0vj4+O1devWRf9BfT1tR1l++eUXjYyM1Pvvv9/rdcaMGaMiUjRwwBeOHDnisx+qyqSlpSmgjz32WI23VVBQoM8995xGRERos2bNNCYmptam4Pjmm2+0Xbt2Gh0drVOmTKm0fn5+vi5atEhTUlKK4uzSpYs+++yz+vPPP1sfRXUelijqhttuu01jY2P10KFDZS4/duyYNmnSRK+55hqf7nfr1q0K6JNPPumT7T311FOlfiBra56f3r17a8+ePb2qu2HDBo2IiND/+7//83kctTmpYOEopIqu/ajM/v37i/oMhg8frvv27av1Sf2ys7O1f//+CujYsWP13//+d6n9b9u2TR955BFNSkpSQOPj4/W2227TFStWlDojtFFPlijqnePHj2vjxo111KhRFda7++67NSoqSvfs2eOzfT/99NMK+Gya9kBOajhu3DgNDw/3quN80KBB2qhRo0pHKlVHbb4H+/bt0zZt2mjHjh31yJEjVV7/q6++0nbt2mlERIQ+99xzfhtJ5Y3c3Fy96667io2yKjnqCtwXoKalpfn92iJLFCaovPvuu+rNhVRr1qxRwKdt3d26ddPzzjvPZ9sL5BTdH330kQL6n//8p8J68+fPV0D/9re/+SWO2p49tXB472233eb1OgUFBfrSSy9pVFSUtm3bVpctW+aX2KqjSZMmZX6H4uPjdevWrbUWhyUKE1SuuOIKbdGihVedkt27d9ezzz7bJ3/5ff/99wros88+W+NtFQrkGcWhQ4c0PDxcx48fX26d3Nxc7dSpk3bo0EGPHz/ut1hqu+mmcLSPN535Bw4c0KuvvloB/fWvf+3TM1RfCOQfG54sUZigsW/fPo2KitJ77rnHq/qFV277YlK6wovFtm/fXuNtFQr0vQjOPfdcvfDCC8td/vzzzytQ42sZgs2xY8f0zDPP1BYtWlTYnLZ69Wpt3769hoeH69NPP13qzpHBIND3ZClkicIEjVdffVWBoiGlldm/f7/GxMT4ZHK4M844Q3v37l3j7ZQUyLub3XvvvRoVFVVm+3V2drYmJCToJZdcEtC2eH/55ptvNCoqSocNG1bq+AoKCvTVV1/V6Ohobd26tX766acBirJygf5jo5AlChM0+vTpo6effnqVfriuu+46TUhIqFFn3vr1633e3xEM5syZo4B+/PHHpZbdfvvtGh4eruvXrw9AZLXj2WefVUAnT55cVHb48OGieakuvfRS3b17dwAj9E4w3ErVEoUJClu2bFFAn3jiiSqtl5GRoYC+9dZb1d73Qw89pCKiO3furPY2gtG+fftURPQvf/lLsfJvvvlGw8LC9He/+12AIqsd+fn5OmDAAI2MjNRWrVqpiGhERIQC+vjjjwdlU1OwqkmiCMMYHym8rWTh7SG91b9/f5KSknj99dertV9VJT09neTkZFq2bFmtbQSrxo0bc9ZZZ/HJJ58UlakqY8aMISEhgUceeSRwwdWCsLAwhg8fzokTJ9i5cyeqSl5eHtHR0SQlJREWZj9htcHeZeMTqsrUqVPp168fiYmJVVo3LCyMm266iUWLFrFt27Yq7/vrr79m48aNjBw5ssrr1gV9+/Zl2bJlRfecfv/991m8eDF/+ctfaNKkSYCj87/nnnuuVNnx48cZP358AKIJTZYojE+sWLGCTZs28Zvf/KZa648ePRoR4c0336zyuunp6YSHh3PllVdWa9/BLjk5maNHj7Jy5UqOHz/OvffeS5cuXbjtttsCHVqt2L59e5XKje9ZojA+MXXqVGJiYrjqqquqtX5iYiIXX3wxb7zxBvn5+V6vV9jsdNFFF9GsWbNq7TvY9e3bF4ClS5fywgsvsGXLFp5//nkiIiICHFntKO8Mtapnrqb6LFGYGjtx4gQzZsxg2LBhxMfHV3s7N998M9u3b2fx4sVer7Ny5Uq2bt3KtddeW+39BrvmzZvTqlUrHn74YR588EFiY2PZvXt3oMOqNRMmTCAuLq5YWVxcHBMmTAhQRKHHEoWpsY8++oi9e/dy/fXX12g7l19+OU2aNKlSp3Z6ejqRkZGMGDGiRvsOZmlpaezevbuoj+Lo0aOkpqaSlpYW4MhqR0pKCpMnT8blciEiuFwuJk+eXOVBE6b6xD1qqu7o0aOHrly5MtBhGA/XXHMNS5YsYefOnURGRtZoW3fffTeTJk3ip59+qrSjtqCggKSkJM466yzmzZtXo/0Gs6SkJDIzM0uVu1yuanX+m9AkIqtUtUd11rUzClMjBw4cYM6cOYwcObLGSQLglltuITc3l7feeqvSusuXLycrK6teNzuBdeaawLNEYWrk7bff5vjx4zVudip09tln061bN6+an2bMmEF0dDTDhw/3yb6DlXXmmkDzOlGISB8Rucl53lxE2vkvLFNXTJ06lY4dO3Luuef6bJs333wza9asYfXq1eXWyc/PZ9asWQwePJhGjRr5bN/ByDpzTaB5lShE5GHgAWCcUxQJTPNXUKZuyMzMZOnSpVx//fWIiM+2e9111xEdHc2//vWvcut8+umn7Nq1q943O4F15prA8/aMYgQwDPgFQFV3Aif5KyhTNxSOuqnuRXblady4MVdccQVpaWkcO3aszDrp6enExcUxZMgQn+47WKWkpLBt2zYKCgrYtm2bJQlTq7xNFLnOpFIKICIN/BeSKSktLa1oXpukpKSgGBZZOGXHhRdeSFJSks+3f8stt7B//37ee++9Usvy8vKYPXs2Q4YMoUED+yoa42/eJoqZIjIJSBCR/wMygH/6LyxTKC0tjdTUVDIzM1FVMjMzg2IM/apVq/juu+98fjZRqH///rhcrjI7tZcsWUJ2dnZINDsZEwy8ShSq+jfgbWA20Al4SFX/4c/AjNv48ePJyckpVpaTkxPwCdGmTp1KVFQUV199tV+2X9FEgTNmzOCkk05i0KBBftm3MaY4bzuz2wGfqur9qnof8JmIJHmx3mUislFENovIg+XUuUZEvhWR9SJS+eD5EBOMY+hPnDjB9OnTGTp0KI0bN/bbfkaPHg3AlClTispyc3N55513GD58OLGxsX7btzHmf7xtepoFFHi8znfKyiXsUVZCAAAaqUlEQVQi4cBEYBDQGRglIp1L1OmAeyRVb1XtAozxMp6QUd5YeVXl1ltvZcOGDbUcESxcuJDs7GyfXTtRHpfLxcCBA3njjTcoKCgo2veBAwes2cmYWuRtoohQ1dzCF87zqErWOQ/YrKpbnPozgJJXRv0fMFFV9zvbDZ2ZzrxU1hj6mJgY+vfvz5QpU+jcuTODBw9m4cKF1NZ0LFOnTqVp06a10vRz8803k5mZWTRRYHp6OgkJCVxyySV+37cxxs3bRJEtIsMKX4jIcGBPJeu0BrI8Xu9wyjx1BDqKyOcislxELvMynpCRkpLCPffcA1A0hv61115j8eLFZGVl8eijj7Jq1SouueQSzj77bN544w2OHz/ut3gOHTrEe++9x7XXXktUVGV/K9Tc5ZdfTuPGjXn99dc5duwY7733HiNGjKiVfRtj3LxNFLcBfxSR7SKShfviu1t9sP8IoAPQDxgF/FNEEkpWEpFUEVkpIiuzs7N9sNu6Jzw8nIMHDxYbQ9+8eXP+/Oc/s337dt544w3A/Re4y+Xiscceo/C98uXw2tmzZ3Ps2DG/jXYqKSYmhpSUFN59912mT5/O4cOHrdnJmNpWlRtsAw2Bhl7W7QUs8Hg9DhhXos6rwE0erxcB51a03e7du1f75uJ11bnnnqu9e/eutF5BQYEuXLhQBw8erIDGxMRo//79NSYmpvAaGAU0Li5Op02bVq1Y+vfvr+3bt9eCgoJqrV8dEyZMKIo9LCxMp0yZUmv7Nqa+AFZqFX7vPR/eJoho4Drgj8BDhY9K1okAtgDtcPdnfA10KVHnMmCK87wZ7qaqphVtN9QSxb59+1RE9OGHH67Set9++62mpqYWSxCeD5fLVeVYtm/friKijzzySJXXra5p06ZpXFyczxKdMaGqJonC26an93F3ROfhnsaj8FHRmUoecBewANgAzFTV9SLyqEd/xwJgr4h8CywB7lfVvV7GFBKWLFmCqjJw4MAqrXf66aczadKkcudgyszM5LbbbuPFF18kIyODnTt3ltsZXth0lZiYiKrW6iR8wXodiTGhxKsbF4nIOlU9oxbiqVSo3bjojjvuYOrUqezbt69a93so76Y30dHRxMXFsX///qKy+Ph4OnfuTOfOnTn99NPp3Lkz33//PX/84x+L/VjHxcXV2qR0YWFhZSYwESkaMmuMqVxNblzkbaKYDPxDVddWZye+FGqJomPHjnTs2LHad3ArnAKkrB/66667jt27d/Ptt98We2zYsIGff/65wu3W1t3V7O5uxvhGTRKFt30U3wK5wEbgG2At8E1127tq8gilPopt27YpoM8//3yNtjNt2jR1uVwqIupyubxq39+zZ49++umn5fZxiEiNYvKW9VEY4xvUoI8iwst8YpPqBMCiRYsAqtw/UVJKSkqVm4maNm1Knz59cLlcZf5FX1t3VyuMe/z48Wzfvp3ExEQmTJhg02wbU4u8nRQwU1UzgaMU/8vS+FFGRgYnn3wyXbp0CVgMwXB3NbsXgzGB5e2kgMNEZBOwFfgE2AZ86Me4Ql5BQQEZGRkMHDjQp3ePqyq7u5oxxtump8eA84EMVT1HRPoDtXNpbohat24d2dnZNW528oXqNF0ZY+oPb6+jOKHu6xvCRCRMVZcA1es9N17JyMgA4KKLLgpwJMaYUOftGcUBEWkILAXSRGQ3lVxwZ2omIyODTp060bZt20CHYowJcd6eUQzH3ZH9e+Aj4AdgqL+CCnW5ubl88sknQdHsZIwxXp1RqKrn2cOUcisan1i+fDk5OTmWKIwxQaHCRCEin6lqHxE5TPHhsAKoqtbepD8hJCMjg7CwMPr16xfoUIwxpuJEoap9nH9Pqp1wDLgTxbnnnktCQqlbcxhjTK2rtI9CRMJF5LvaCMbAwYMHWbFihTU7GWOCRqWJQlXzgY0iUjtzNoS4Tz75hPz8fEsUxpig4e3w2MbAehFZgcewWFUdVv4qpjoyMjKIjY2lV69egQ7FGGMA7xPFn/0ahSmSkZFB3759iY6ODnQoxhgDeD889hN/B2Lgxx9/ZMOGDdx8882BDsUYY4p4Oyng+SLylYgcEZFcEckXkUP+Di7U+GpacWOM8SVvr8x+CRgFbAJigd8CE/0VVKjKyMigWbNmnHXWWYEOxRhjinibKFDVzUC4quar6hvAZf4LK/SoKhkZGVx00UWEhXn9sRhjjN9525mdIyJRwBoReQb4iSokGVO5DRs28NNPP1mzkzEm6Hj7Y3+9U/cu3MNj2wJX+iuoUFQ4rfjFF18c4EiMMaY4b88ougMfqOoh4C9+jCdkZWRk0L59e1wuV6BDMcaYYrw9oxgKfC8iU0VkiIh4m2CMF06cOMHHH39szU7GmKDkVaJQ1ZuA9sAs3KOffhCR1/wZWCj56quvOHz4sCUKY0xQ8vrMQFVPiMiHuKcbjwUuxz1M1tRQRkYGIkL//v0DHYoxxpTi7QV3g0TkTdzXUVwJvAb8yo9xhZSMjAy6d+9OkyZNAh2KMcaU4m0fxQ3Ae0AnVR2tqvNVNc+PcYWMI0eOsGzZMmt2MsYELW/nehrl70BC1dKlS8nLy7NEYYwJWhWeUYjIZ86/h0XkUMl/ayfE+i0jI4OYmBh69+4d6FCMMaZMdivUAMvIyKBPnz7ExMQEOhRjjCmT16OeRKQx7iuyi9ZR1f/6I6hQsWvXLtauXUtKSkqgQzHGmHJ5lShE5DFgNLAFKHCKFRjgn7BCw+LFiwGbVtwYE9y8HfV0DXCqqiaran/nUWmSEJHLRGSjiGwWkQcrqHeliKiI9PA28PogIyODJk2a0LVr10CHYowx5fI2UawDEqqyYREJx33PikFAZ2CUiHQuo95JwD3Al1XZfl1XOK34gAEDCA8PD3Q4xhhTLm8TxZPAahFZICJzCh+VrHMesFlVt6hqLjADGF5GvceAp4FjXkddD2zatImsrCxrdjLGBD1vO7On4P4xX8v/+igq0xrI8ni9A+jpWUFEugFtVfUDEbm/vA2JSCqQCpCYmOjl7oNb4bTiliiMMcHO6xsXqeqLvtyxiIQBz+HuJK+Qqk4GJgP06NFDfRlHoGRkZJCUlMQpp5wS6FCMMaZC3iaKT0XkSWAOcLywsJLhsT/iHk5bqI1TVugk4AzgYxEB99xRc0RkmKqu9DKuOik/P5/Fixdz9dVX4xy7McYELW8TxTnOv+d7lFU2PPYroIOItMOdIEYC1xWtrHoQaFb4WkQ+Bu6r70kCYNWqVRw8eNCanYwxdYK3cz1Vef5rVc0TkbuABUA48C9VXS8ijwIrVbWyzvB6q7B/YsAAuwzFGBP8vL3g7mTgCaCVqg5yhrn2UtXXK1pPVecD80uUPVRO3X5eRVwPZGRk0LVrV5o3bx7oUIwxplLeDo99E/eZQSvn9ffAGH8EVN/l5OTw+eefW7OTMabO8DZRNFPVmThDY517UeT7Lap67LPPPiM3N9cShTGmzvA2UfwiIk1xd2AjIucDB/0WVT2WkZFBVFQUffr0CXQoxhjjFW9HPY3FPTT2VBH5HGgOXOW3qOqxjIwMLrjgAho0aBDoUIwxxiuV3bjoXBH5lXO9RDLwR9zXUfwH95XWpgr27NnD6tWrrdnJGFOnVNb0NAnIdZ5fAIzHPdHffpwrpY33bFpxY0xdVFnTU7iq7nOeXwtMVtXZwGwRWePf0OqfjIwM4uPj6d69e6BDMcYYr1V2RhEuIoXJ5CJgsccyr++OF+rS0tJISkrin//8J7m5uaSnpwc6JGOM8VplP/bTgU9EZA9wFPgUQETaY6OevJKWlkZqaio5OTkAHD16lNTUVAC7Baoxpk4Q1YonY3WGwrYE/qOqvzhlHYGGgbhndo8ePXTlyrozHVRSUhKZmZmlyl0uF9u2bav9gIwxIUlEVqlqte4iWmnzkaouL6Ps++rsLBRt3769SuXGGBNsvL3gzlSDqtKwYcMyl9WXGzAZY+o/SxR+oqqMGTOGw4cPExFR/MQtLi6OCRMmBCgyY4ypGksUflCYJF588UXGjh3Lm2++icvlQkRwuVxMnjzZOrKNMXWGDXH1MVXlnnvu4R//+Af33nsvf/3rXxERSwzGmDrLzih8SFW5++67+cc//sF9991XlCSMMaYuszMKH1FVfve73zFx4kTuu+8+nnnmGUsSxph6wc4ofMAzSdx///2WJIwx9YolihpSVe666y4mTpzIH/7wB55++mlLEsaYesUSRQ0UFBRw55138vLLL/OHP/yBp556ypKEMabesURRTQUFBdx111288sorPPDAA5YkjDH1liWKaig8k3jllVd48MEHefLJJy1JGGPqLRv1VEUFBQXccccdTJo0iQcffJAnnnjCkoQxpl6zMwovFN5PIiwsjPj4eCZNmsS4ceMsSRhjQoKdUVSi5P0kjhw5QkREBF26dLEkYYwJCXZGUYnx48cXJYlCeXl5jB8/PkARGWNM7bJEUQm7n4QxJtRZoqhEefeNsPtJGGNChSWKSkyYMIHIyMhiZXY/CWNMKLFEUYmUlBTatGlDVFSU3U/CGBOSbNRTJbZs2cLWrVt55plnuP/++wMdjjHG1Do7o6jErFmzALjmmmsCHIkxxgSGXxOFiFwmIhtFZLOIPFjG8rEi8q2IfCMii0TE5c94qiM9PZ2ePXvicgVdaMYYUyv8lihEJByYCAwCOgOjRKRziWqrgR6qehbwNvCMv+Kpjk2bNrF69WquvfbaQIdijDEB488zivOAzaq6RVVzgRnAcM8KqrpEVQuvZlsOtPFjPFWWnp4OwNVXXx3gSIwxJnD8mShaA1ker3c4ZeW5BfiwrAUikioiK0VkZXZ2tg9DrNjMmTPp06cPbdoEVf4yxphaFRSd2SLyG6AH8NeylqvqZFXtoao9mjdvXisxbdiwgbVr11ontjEm5PlzeOyPQFuP122csmJEZCAwHkhW1eN+jKdK0tPTERGuuuqqQIdijDEB5c8ziq+ADiLSTkSigJHAHM8KInIOMAkYpqq7/RhLlagq6enpJCcn07Jly0CHY4wxAeW3RKGqecBdwAJgAzBTVdeLyKMiMsyp9legITBLRNaIyJxyNler1q1bx3fffWejnYwxBj9fma2q84H5Jcoe8ng+0J/7r6709HTCwsK44oorAh2KMcYEXFB0ZgeTwmanAQMG0KJFi0CHY4wxAWeJooTVq1ezefNma3YyxhiHJYoSZs6cSUREBCNGjAh0KMYYExQsUXgobHYaOHAgTZs2DXQ4xhgTFCxRePjqq6/Ytm2bNTsZY4wHSxQe0tPTiYqK4vLLLw90KMYYEzQsUTgKCgqYNWsWl156KQkJCYEOxxhjgoYlCsfy5cvJysqyuZ2MMaYESxSO9PR0oqOjGTZsWOWVjTEmhFiiAPLz85k1axaDBw+mUaNGgQ7HGGOCiiUK4PPPP+enn36y0U7GGFMGSxS4m51iY2P59a9/HehQjDEm6IR8osjLy+Ptt99myJAhNGzYMNDhGGNM0An5RPHJJ5+we/dua3YyxphyhHyimDlzJg0aNGDw4MGBDsUYY4JSSCeKEydOMHv2bIYNG0ZsbGygwzHGmKAU0oli8eLF7N2715qdjDGmAiGdKNLT02nUqBGXXXZZoEMxxpigFbKJIjc3l3fffZfLL7+c6OjoQIdjjDFBK2QTxcKFCzlw4IDN7WSMMZUI2USRnp5O48aNufjiiwMdijHGBLWQTBTHjh3jvffeY8SIEURFRQU6HGOMCWohmSgWLFjA4cOHbbSTMcZ4ISQTRXp6Ok2bNqV///6BDsUYY4JeyCWKnJwc5syZw5VXXklkZGSgwzHGmKAXconiww8/5JdffrFmJ2OM8VLIJYr09HRatGhBcnJyoEMxxpg6IaQSxZEjR5g3bx5XXXUV4eHhgQ7HGGPqhJBKFPPmzePo0aPW7GSMMVUQEokiLS2NpKQkRo0aRXh4ONu3bw90SMYYU2dEBDoAf0tLSyM1NZWcnBwA8vPzufXWWxERUlJSAhydMcYEv3p/RjF+/PiiJFEoJyeH8ePHBygiY4ypW+p9oiivmcman4wxxjt+TRQicpmIbBSRzSLyYBnLo0Uk3Vn+pYgk+TqGxMTEKpUbY4wpzm+JQkTCgYnAIKAzMEpEOpeodguwX1XbA88DT/s6jgkTJhAXF1esLC4ujgkTJvh6V8YYUy/584ziPGCzqm5R1VxgBjC8RJ3hwBTn+dvARSIivgwiJSWFyZMn43K5EBFcLheTJ0+2jmxjjPGSP0c9tQayPF7vAHqWV0dV80TkINAU2ONZSURSgVSoXpNRSkqKJQZjjKmmOtGZraqTVbWHqvZo3rx5oMMxxpiQ4s9E8SPQ1uN1G6eszDoiEgHEA3v9GJMxxpgq8mei+AroICLtRCQKGAnMKVFnDnCj8/wqYLGqqh9jMsYYU0V+66Nw+hzuAhYA4cC/VHW9iDwKrFTVOcDrwFQR2Qzsw51MjDHGBBG/TuGhqvOB+SXKHvJ4fgy42p8xGGOMqRmpay09IpINZJazuBklRkyFGDv+0D5+sPfAjr/843eparVGA9W5RFEREVmpqj0CHUeg2PGH9vGDvQd2/P45/joxPNYYY0zgWKIwxhhTofqWKCYHOoAAs+M3of4e2PH7Qb3qozDGGON79e2MwhhjjI9ZojDGGFOhepMoKrtJUl0lIm1FZImIfCsi60XkHqe8iYgsFJFNzr+NnXIRkRed9+EbEenmsa0bnfqbROTG8vYZbEQkXERWi8g853U750ZXm50bX0U55eXeCEtExjnlG0Xk0sAcSfWISIKIvC0i34nIBhHpFWKf/++d7/46EZkuIjH1+TsgIv8Skd0iss6jzGeft4h0F5G1zjovinhxawdVrfMP3FOE/ACcAkQBXwOdAx2Xj46tJdDNeX4S8D3uG0E9AzzolD8IPO08Hwx8CAhwPvClU94E2OL829h53jjQx+flezAWeAuY57yeCYx0nr8K3O48vwN41Xk+Ekh3nnd2vhPRQDvnuxIe6OOqwvFPAX7rPI8CEkLl88d9K4KtQKzHZz+6Pn8HgL5AN2CdR5nPPm9ghVNXnHUHVRpToN8UH72xvYAFHq/HAeMCHZefjvV94GJgI9DSKWsJbHSeTwJGedTf6CwfBUzyKC9WL1gfuGcdXgQMAOY5X+49QETJzx73vGK9nOcRTj0p+X3wrBfsD9wzKm/FGXhS8nMNgc+/8J41TZzPdB5waX3/DgBJJRKFTz5vZ9l3HuXF6pX3qC9NT2XdJKl1gGLxG+c0+hzgS+BkVf3JWbQLONl5Xt57UVffoxeAPwAFzuumwAFVzXNeex5HsRthAYU3wqqrxw7uv36zgTec5rfXRKQBIfL5q+qPwN+A7cBPuD/TVYTWdwB893m3dp6XLK9QfUkU9Z6INARmA2NU9ZDnMnX/aVDvxjmLyBBgt6quCnQsARSBuxniFVU9B/gFd9NDkfr6+QM4bfHDcSfMVkAD4LKABhVggfi860ui8OYmSXWWiETiThJpqvqOU/yziLR0lrcEdjvl5b0XdfE96g0ME5FtuO+5PgD4O5Ag7htdQfHjKO9GWHXx2AvtAHao6pfO67dxJ45Q+PwBBgJbVTVbVU8A7+D+XoTSdwB893n/6DwvWV6h+pIovLlJUp3kjEh4Hdigqs95LPK86dONuPsuCstvcEZDnA8cdE5ZFwCXiEhj56+0S5yyoKWq41S1jaom4f5MF6tqCrAE942uoPSxl3UjrDnASGdETDugA+4OvaCnqruALBHp5BRdBHxLCHz+ju3A+SIS5/xfKDz+kPkOOHzyeTvLDonI+c77eYPHtsoX6E4bH3b+DMY9IugHYHyg4/HhcfXBfZr5DbDGeQzG3e66CNgEZABNnPoCTHTeh7VAD49t3Qxsdh43BfrYqvg+9ON/o55Owf2ffDMwC4h2ymOc15ud5ad4rD/eeU824sUoj2B6AF2Blc534D3co1hC5vMH/gJ8B6wDpuIeuVRvvwPAdNz9MSdwn1He4svPG+jhvJc/AC9RYqBEWQ+bwsMYY0yF6kvTkzHGGD+xRGGMMaZCliiMMcZUyBKFMcaYClmiMMYYUyFLFKZeEZEjftjmNme2zbXinsX3cRGJcZa1EpG3K1g3QUTu8HVMxtQmGx5r6hUROaKqDX28zW24x6fvcaZSmQycUNVKp+p25ueap6pn+DImY2qTnVGYek9Ehjr3JlgtIhkicrJT3tyZ23+9M9lepog0q2hbqnoEuA243LlHQFLhfQNEpIuIrBCRNc69AToATwGnOmV/FZGGIrJIRP7rnKEMd9ZNEve9Jv7pxPMfEYl1lrV34v7aWe9Up/x+EfnK2ddf/PcOmlBnicKEgs+A89U9qd4M3LPRAjyMe4qHLrjnUEr0ZmPqnpRxK+5pIDzdBvxdVbvivvp1B+4J/H5Q1a6qej9wDBihqt2A/sCzHjeO6QBMdOI5AFzplKc55WcDFwA/icglTv3zcF+53V1E+nr9jhhTBRGVVzGmzmsDpDuTqUXh/pEH9/QoIwBU9SMR2V+FbZZ1V7BlwHgRaQO8o6qbyrh5mABPOD/qBbineC6cMnqrqq5xnq8CkkTkJKC1qr7rxHkMwEkUlwCrnfoNcSeOpVU4BmO8YmcUJhT8A3hJVc8EbsU9H1C1OT/eSbjnFiuiqm8Bw4CjwHwRGVDG6ilAc6C7c+bxs0c8xz3q5VPxH3ICPOmcqXRV1faq+np1jseYyliiMKEgnv9NpezZAf05cA0U/YXeuLINOZ3ZLwPvqer+EstOAbao6ou4Z+Q8CziM+xa2nrHsVtUTItIfcFW0P1U9DOwQkcudfUSLSBzu2UFvduJBRFqLSIvK4jemOqzpydQ3cSLieQev54BHgFlO09Ji3DfBAfespNNF5HrczUa7cP+wl2WJ05cQBrwLPFZGnWuA60XkhLOtJ1R1n4h87nR4fwg8DcwVkbW4Z4T9zotjuh6YJCKP4p5R9GpV/Y+InA4sc5q3jgC/4X/3KTDGZ2x4rAlZIhIN5Ktqnoj0wn0Xua6BjsuYYGNnFCaUJQIzRSQMyAX+L8DxGBOU7IzCGGNMhawz2xhjTIUsURhjjKmQJQpjjDEVskRhjDGmQpYojDHGVOj/AZak9pVh36qEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "geoplot.semivariogram(P, lags, tolerance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can pass a model to this function using the optional `model` argument and see it plotted in red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4FNX6wPHvm0pCC02kJihFwIKCIopEEVFQQK6KYK6K5UbFhmC92AF7Qa+ocPWqP4gQ7IAgGlCxgBRBAZEiEEBUAqEHSHt/f8wkpmeT7GY3yft5nn2yO3Nm5t2SfXfOOXOOqCrGGGNMcYL8HYAxxpjAZonCGGNMiSxRGGOMKZElCmOMMSWyRGGMMaZEliiMMcaUyBKFqdFE5BwRWVfV9m1MZRK7jsIEKhHpCTwDdAaygLXASFVd6tfAjKlhQvwdgDFFEZF6wGzgFmAGEAacAxz1Z1yeEpEQVc2savs2pihW9WQCVXsAVZ2mqlmqelhVP1fVnwFE5HoRWSsie0RknohE52woIioiI0Rkg4gcEJGxInK8iHwvIvtFZIaIhLllzxWR7e79+0Tk/bxBiMhLIvKye/8695gHRGSTiNyUp9y5IrLd3cefwFt59+2W6SgiX4nIXhFZIyID86xrJCKz3PiWisg4Efm2wHO6VUQ2ABvyxLbN3Wa5iJyTp/yjIvKeiEx1410lIu1F5AER2elu19c7b5Wp7ixRmEC1HsgSkXdEpJ+INMhZISKDgH8D/wCaAN8A0wpsfyHQFTgTuBeYDPwTaAWcCAwr4pjTgf4iUtc9TjAwBHjXXb8TuASoB1wHvCgip+XZ/ligIRANxOfdsYiEArOAz4FjgNuBBBHp4BaZCBxy93GteyvoUqA70Ml9vBTo4h7zXeA9EamVp/wAYArQAFgBzMP5n28BPA5MKuIYxhSmqnazW0DegI7A28B2IBOYCTQF5gI35CkXBKQB0e5jBc7Os345cF+ex88DE9z75wLb86z7FrjGvX8B8FsJ8X0M3JlnP+lArTzrc/eNU232JxCUZ/004FEgGMgAOuRZNw74Ns9jBXqX8nrtAU5x7z8KfJFn3QDgIBDsPq7r7jPK3++z3QL/ZmcUJmCp6lpVHa6qLXHOApoDE3B+sb/kVuHsBVIBwfmlnOOvPPcPF/G4TjGHfZe/zzau4u+zCdwzm8Uikuoetz/QOM+2Kap6pJj9Nge2qWp2nmXJbsxNcNoLt+VZl/d+kctE5G63KmyfG0/9AvEUfM67VDUrz2Mo/nUwJpclClMlqOqvOGcXJ+J8Yd6kqlF5bhGq+r0XDvUecK6ItAQG4yYKEQkHPgCeA5qqahQwBydB5YZZwn53AK1EJO//XGvgdyAF54ypZZ51rYrYR+7+3faIe3Gqxhq48ewrEI8xXmGJwgQkETlBREa7X9iISCucX/qLgdeBB0Sks7uuvohc4Y3jqmoK8BXwFrBZVde6q8KAcNwvdRHpB5SlMfgHnOqxe0UkVETOxakOmu7+yv8QeFREIkXkBOCaUvZXFye5pAAhIvIwTtuJMV5nicIEqgM4Dbc/iMghnASxGhitqh8BTwPTRWS/u7yfF4/9LtCHPNVOqnoAuAOnq+4enGqpmZ7uUFXTcRJDP2AX8CpOW8ivbpHbcKqO/sRpgJ5GyV2B5wGf4TT6JwNHKLq6ypgKswvujAlAIvI0cKyqFtX7yZhKZWcUxgQAt6rtZHGcAdwAfOTvuIwBuzLbmEBRF6e6qTlOb6XngU/8GpExLqt6MsYYUyKrejLGGFOiKlf11LhxY42JifF3GMYYU6UsX758l6o2Kc+2VS5RxMTEsGzZMn+HYYwxVYqIJJd3W59VPYnI/9xRKlcXs15E5GUR2SgiPxcYXM0YY0yA8GUbxdvARSWs7we0c2/xwGs+jMUYY0w5+SxRqOpCnMHaijMI+D91LAaiRKSZr+IxxhhTPv7s9dSC/EMObCf/6J+5RCReRJaJyLKUlJRKCc4YY4yjSnSPVdXJqtpNVbs1aVKuRntjjDHl5M9E8Tv5h1Ju6S4zplpJSEggJiaGoKAgYmJiSEhI8HdIxpSJPxPFTOAat/fTmcA+Vf3Dj/EY43UJCQnEx8eTnJyMqpKcnEx8fLwlC1Ol+GwIDxGZhjMVZGOcsWseAUIBVPV1ERHgFZyeUWnAdapa6gUS3bp1U7uOwlQVMTExJCcX7r4eHR3Nli1bKj8gU2OJyHJV7VaebX12wZ2qFjV5fd71Ctzqq+MbEwi2bt1apuXGBKIq0ZhtTFXVunXrMi03JhBZojDGh8aPH09kZGS+ZZGRkYwfP95PERlTdpYojPGhuLg4br01fw3riy++SFxcnJ8iMqbsLFEY42MHDx4kMjKSuXPnAtC8eXM/R2RM2ViiMMaHVJVZs2bRt29fYmNjCQsL4+uvv/Z3WMaUiSUKY3xo5cqVbN++nQEDBhAREcEZZ5zBwoUL/R2WMWViicIYH5o1axYiwsUXXwxAbGwsy5cv58CBA36OzBjPWaIwxodmzZrFGWecQdOmTQHo1asXWVlZLFq0yM+RGeM5SxTG+MiOHTtYtmwZAwcOzF121llnERwcbO0UpkqxRGGMj3z66acADBgwIHdZnTp16Nq1q7VTmCrFEoUxPjJz5kyio6M58cQT8y2PjY1lyZIlHD582E+RGVM2liiM8YG0tDSSkpIYMGAAzviXf4uNjSU9PZ3Fixf7KTpjysYShTE+MH/+fI4cOZKvfSLH2WefjYhY9ZOpMixRGOMDs2bNom7dusTGxhZaFxUVRZcuXaxB21QZliiM8bLs7Gxmz57NhRdeSFhYWJFlevXqxaJFi0hPT6/k6IwpO0sUxnjZjz/+yB9//JGvt1NBsbGxHDlyhKVLl1ZiZMaUjyUKY7xs5syZBAUF0b9//2LLnHPOOQDWTmGqBEsUxnjZrFmzOOuss2jcuHGxZRo3bkznzp2tncJUCZYojPGibdu2sXLlyhKrnXL06tWL7777jszMzEqIzJjys0RhjBfNnj0bwKNEERsby8GDB1mxYoWvwzKmQixRGONFs2bN4vjjj+eEE04otWyvXr0ArPrJBDxLFMZ4ycGDB5k/f36RV2MXpVmzZrRr184atE3As0RhjJd88cUXpKenF3k1dnF69erFN998Q1ZWlg8jM6ZiLFEY4yWzZs2ifv369OzZ0+NtYmNj2bt3L6tXr/ZhZMZUjCUKY7wgOzubTz/9lH79+hEaGurxdjlDfFg7hQlkliiM8YIlS5awc+dOj3o75dW6dWuio6OtncIENEsUxnjBzJkzCQ4Opl+/fmXeNjY2loULF6KqPojMmIqzRGGMF8yaNYtzzjmHBg0alHnbXr16kZKSwq+//uqDyIypOEsUxlTQli1bWL16dZmrnXJYO4UJdJYoyuGzzz6jQ4cOtG3blqeeeqrQ+rvuuosuXbrQpUsX2rdvT1RUFABffvll7vIuXbpQq1YtPv74Y8CZ6Oa0006jS5cu9OzZk40bN1bqczLlN2vWLMCzq7GLcvzxx9OsWTNrpzCBS1Wr1K1r167qT5mZmXrcccfpb7/9pkePHtWTTz5Z16xZU2z5l19+Wa+77rpCy3fv3q0NGjTQQ4cOqapqu3bt9JdfflFV1YkTJ+q1117rk/iN911wwQXaoUOHCu1j6NCh2rx5c83OzvZSVMbkByzTcn7v2hlFGS1ZsoS2bdty3HHHERYWxtChQ/nkk0+KLT9t2jSGDRtWaPn7779Pv379iIyMBEBE2L9/PwD79u2jefPmvnkCNVBCQgIxMTEEBQURExNDQkKC1/a9f/9+vvrqq3KfTeTo1asXO3bs4LfffvNSZMZ4T4gvdy4iFwEvAcHAG6r6VIH1rYF3gCi3zP2qOseXMVXU77//TqtWrXIft2zZkh9++KHIssnJyWzevJnevXsXWjd9+nRGjRqV+/iNN96gf//+REREUK9ePRYvXuz94GughIQE4uPjSUtLA5z3JD4+HoC4uLgK73/evHlkZGSU6WrsouS0UyxcuJC2bdtWOC5jvMlnZxQiEgxMBPoBnYBhItKpQLEHgRmqeiowFHjVV/H4w/Tp07n88ssJDg7Ot/yPP/5g1apVXHjhhbnLXnzxRebMmcP27du57rrr8iURU35jxozJTRI50tLSGDNmjFf2P2vWLBo2bEiPHj0qtJ+OHTvSuHFja9A2AcmXVU9nABtVdZOqpgPTgUEFyihQz71fH9jhw3i8okWLFmzbti338fbt22nRokWRZadPn15ktdOMGTMYPHhw7hW8KSkp/PTTT3Tv3h2AK6+8ku+//94H0dc8W7duLdPyssjKymLOnDn079+fkJCKnZyLCL169fJZg7Yvq9+qwvFNBZW3caO0G3A5TnVTzuOrgVcKlGkGrAK2A3uArsXsKx5YBixr3bq1l5t4yiYjI0PbtGmjmzZtym3MXr16daFya9eu1ejo6CIbJ7t3764LFizIt89GjRrpunXrVFX1jTfe0H/84x++exI1SHR0tOL8IMl3a968eYX3/c033yigiYmJXohUdcKECQpocnKyV/aXY+rUqRoZGZnv+UdGRurUqVO9epxAPb5xUIHGbH8nilHAaPd+D+AXIKik/fq715Oq6qeffqrt2rXT4447TseNG6eqqg899JB+8sknuWUeeeQRve+++wptu3nzZm3evLlmZWXlW/7hhx/qiSeeqCeffLLGxsbqb7/95tsnUYmmTp2q0dHRKiIaHR1dqV8QU6dO1bCwsEKJokmTJpqSklKhfd9zzz0aEhKie/fu9UqsK1asUECnTJnilf3lKC5ZRkdHe/U4gXp84wjURNEDmJfn8QPAAwXKrAFa5Xm8CTimpP0GQqIwnguEX5OnnnqqBgcH5yaqRx99VGvVqqVnn322HjlypNz7PeGEE7RPnz5eizMzM1OjoqL0xhtv9No+VVVFpMgvahHx6nGKkpmZWeSxK+v43uLPHzveEqiJIsT94m8DhAE/AZ0LlJkLDHfvd8Rpo5CS9muJomrx96/J1NRUDQ0N1bvuuivf8sTERAU0Li6uXNcubNiwQQF96aWXvBWqqqpecskl2r59e6/us2nTpkW+B76uxj1w4IAOGDCg2ERRVc4oAuHHjjdUJFGIs71viEh/YAJO19f/qep4EXncDXim2wvqv0Ad9w24V1U/L2mfdevW1a5du+ZbNmTIEEaMGEFaWhr9+/cvtM3w4cMZPnw4u3bt4vLLLy+0/pZbbuHKK69k27ZtXH311YXWjx49mgEDBrBu3TpuuummQusffPBB+vTpw8qVKxk5cmSh9U888QRnnXUW33//Pf/+978LrZ8wYQJdunQhKSmJcePGFVo/adIkOnTowKxZs3j++ecLrZ8yZQqtWrUiMTGR1157rdD6999/n8aNG/P222/z9ttvF1o/Z84cIiMjefXVV5kxY0ah9V999RUAzz33XO6c0DkiIiKYO3cuAGPHjmX+/Pn51pfUiyc2NpaWLVsydepUAEaOHMnKlSvzlWnfvj2TJ08GIP5f/2L9r79CZiZkZUFmJl1at2bCVVfBvn3887XX2L57d+46MjOprcqc1FR+aNmSp1NS2J2dnbvvrVlZbM7O5rHatXm4Th367dnD4QL/D5eEh3N37doAnJuamrt8e1YWv+XZNk2V/nv2FHqOwyMiGB4Rwa7sbC7fu7fQ+lsiI7myVi22ZWVx9b59bMvKYlN2Nj1CQggTYXTt2gwID2ddZiY3udfZ5PVg7dr0CQ9nZUYGIw8cKLR+XJ06XL13L1uK+D+vC8xr0IAeYWEkHT3KuEOHCpWZVK8eHUJCmHX0KM8XsX5K/fq0Cg4m8cgRXsvTuyxdlVWZmRwCrgkPZ9rRo6Tn2U6ADsHBLG/UiEgRXk1LY8aRI4X2/1XDhgA8d+gQs48ezbcuQoS57thaYw8eZH56er71jYKC+MAdFeGBAwdYlJGRb33L4GCm1q8PwMgDB1hZYH37kBAm16tHTEoKyXk+NznqAvubNgXgn/v2sb3A5FM9QkN5sm5dAC7buzffZw/g/LAwHqpTB8Dzz95xx8GxxwJl/977+uuvl6tqt0IFPeDT6yjUuSZiToFlD+e5/wtwti9jMP4VGRlZqHsqQHh4uOc7WboUJkyA6dOh4D/sihVQ8ILH4GAICYGQENZmZtImIoLTe/aExYshz5dNa6B+ejqP7N5N2+7dYcMGJ8HkFR0Np5zi3J85M3fx7t27iczOpvGpp0LnzpCRAW7CzKdDB+d2+DB88UXh9Z06Qdu2cPAgLFhAVEYG7NrF3jp1OCYiAk4+GWJiYO9eKKpH1GmnQcuWsGsXFNFTbm7TpmzZs4fBMTHM2bqVo9nZhAcFUT8sjJ1HjnCDKl9ccAHs2QM//lh4/716QVQUbNkCP/9ceH3v3lCnDmzcCL/8AsDBjAxWpaaSKUJC794Ma9uW8IULeWfdOo5mZxMEZAORDRrAJZdAaCisWQNFXWyYc33KTz9BcnL+dSEhkPMFuXw5/P57/vW1akHfvs79H36Av/7Kv752bTj/fOf+d9/B7t35Vmv9+nwSHU3yvHmF4wIOAFx6qfNg/nwomEibNgW3JyOffw4FE2GLFpDzo3fOHM8+exERRcbic+U9FfHXzaqeqpaiGpNr1apV+ml7Robqe++pnnWWU0Nat67qLbeoPvOM6qRJqtOnq86dq/rdd6qrV6tu26a6f79qnk4CKSkpGhwcXGSnghxHjx7V2NhYDQsL02+//daj55SamqrBwcH6wAMPeFS+LDIyMrR27do6YsSICu9r37592rRpU+3Ro0eR1WszZ87UOnXqaPPmzXX58uUVPp6q6pw5c7ROnTraokULXblyZZFldu/erS1atNATTjghdwibQHL48GGdPHmydujQQQENDg72S9WdtxGIbRS+ulmiqHpOO+00DQoKym1ULbEBeM8e1eeeU42Odj6exx2n+tJLqvv2lfm4kyZNUkB//PHHEsvt2rVL27Vrp40bN/aot9m7776rgH7//fdljskTffv21c6dO1d4P/fdd58CumTJkmLL/PTTT9q6dWuNiIjQDz74oELHe/XVVzUoKEhPPfVU3b59e4llv/jiCwX01ltvrdAxvWn37t06bty43Dad0047TadNm6b/93//V6iNAtAuXbrovnJ8Lv3FEoUJWLt27dLQ0FAdNWqUqqrGxcVp/fr1NS0tLX/BDRtUb79dtXZt52MZG6v68ceqmZnlPnbv3r21Xbt2HjVWr1+/Xhs2bKgnnHCC7tmzp8Syw4YN0yZNmmhmBWIryfjx4xWoUPfdjRs3alhYmF5zzTWllv3zzz/1zDPPVEDHjx9f5sb9zMxMHTVqlAJ6ySWX6IEDBzzabuTIkQronDlzynQ8b9u8ebPecccdWrt2bQW0X79+On/+/HyvQ95eT61bt9Zrr71Wg4ODtWPHjrp+/Xo/Ru85SxQmYE2cOFEBXbFihaqqzp8/XwFNSEhQzc5WXbBAdeBAVRHV0FDVa65R9UI1yJ9//qlBQUH64IMPerzNV199paGhoXr++edrenp6kWXS09M1KipKhw8fXuEYi5NzId+HH35Y7n0MHjxYa9eurb///rtH5Q8fPqxXXXWVAnr11Vd73G344MGDOmjQIAX0jjvuKFPyPHz4sHbu3FmPPfbYCl/TUpqiurcuW7ZMhw4dqsHBwRoaGqrXXnutrlq1yuN9LliwQBs1aqRRUVE6b948H0bvne65lihMwOrRo4eeeOKJub/OsrKytE1MjPbu2FH1lFOcj2DjxqoPPaS6Y4fXjvvKK68oUKZ/fFXVt99+WwG98cYbi/xl/eWXXypQ4Wqakhw5ckRr1aqlI0eOLNf2Ocl4/PjxZdouOztbx44dq4CeddZZ+tdff5VYfseOHdq1a1cNCgoqdzfhlStXalhYmF566aU+G2K9qO6tQUFBCmi9evX0nnvuKbWqrDibNm3Sk046SYOCgvT555/3yXPwVvdcSxQmIOVca/D00087C/76S/XRR/Vx9xR/U7t2qm+8oVqwGsoLevbsWe56/jFjxiigzz77bKF1o0aN0rCwMI+rV8rr3HPP1VNPPbXM22VkZOhJJ52kMTExevjw4XIde8aMGVqrVi2NiYkpNtH+/PPP2qpVK61du7bOmjWrXMfJ8cwzzyigb775ZoX2U5ziruWJioryylX1Bw4c0MGDByug11xzTblf9+K0bt3aK9ehWKIwAemRRx5REdFtW7eqvvKKani4KujW885TEdGHH3rIJ8fdtm2bAvr444+Xa/usrCwdMmSIiki+6p/s7Gxt27atXnTRRd4KtVg5r11Zv8hee+01BfS9996r0PGXLFmizZo107p16+qnn36ar+rjmGOO0fDwcG3evHmpHQU8kZmZqeeee67WqVNHN27cWOH9FVQZV6ZnZWXpo48+qoB2797d4yq/kqxZs0ZHjx7ttSvbLVGYgJOdna3HH3+89j73XNX4eOejdvHFqr/+qqqqF154obZq1conDcIvvPCCArmDLJZHWlqadu/eXSMiInTZsmWq6gz0COjEiRO9FWqxFixYoIDOnj3b421SU1O1UaNGGhsb65UqkG3btumpp56qgIaGhhb6knr55ZcrfIwcycnJWr9+fe3Ro4dmZGR4bb+LFi0qtnurL64M/+CDD7R27drarFkzXbx4cZm337t3r77++uvavXt3BTQkJEQjIiLsjKKsN0sUVcN3332ngL7Vvr3zMXvggXzXOOQMoeGLRsDu3btrly5dKryfP//8U6Ojo7V+/fraokWL3H9Qbw/bUZRDhw5paGio3nPPPR5vM3LkSBWR3I4D3nDw4EGvfVGVJiEhQQEdO3ZshfeVnZ2tL7zwgoaEhGjjxo21Vq1aFa7j99TPP/+sbdq00fDwcH3nnXdKLZ+VlaXz58/XuLi43Dg7d+6szz//vP7111/WRlGemyWKquHmK67QCBHdHx6umpBQaP2RI0e0YcOGOmTIEK8ed/PmzQrok08+6ZX9PfXUU4W+ICtrnJ+zzz5bu3fv7lHZtWvXakhIiP7rX//yehyVOahgTi+kkq79KM2ePXty2wwGDRqkqamplT6oX0pKip533nkK6KhRo/T//u//Ch1/y5Yt+uijj2pMTIwCWr9+fb355pt1yZIlhc4IrdeTJYpq52hiojYAHRYRoVrCP/wdd9yhYWFhumvXLq8d++mnn1bAa8O0+3NQwwceeECDg4M9ajjv16+f1qtXr9SeSuVRma9BamqqtmzZUtu3b68HDx4s8/ZLly7VNm3aaEhIiL7wwgs+60nlifT0dL3tttvy9bIq2OsKnAtQExISCl9b5GWWKExgyM5WHTtWP3L/AeaUMq/CypUrFfBqXfdpp52mZ5xxhtf2588huj/77DMF9PPPPy+x3Jw5cxTQ5557zidxVPboqTnde2+++WaPt8nOztZXXnlFw8LCtFWrVrpo0SKfxFYeDRs2LPIzVL9+fd28eXOlxWGJwvjfoUOqQ4aogv6jVSs95phjPGqU7Nq1q55yyile+eW3fv16BfT555+v8L5y+POMYv/+/RocHKxjxowptkx6erp26NBB27Vrp0ePHvVZLJVddZPT28eTxvy9e/fqFVdcoYBefPHFXj1D9QZ//tjIyxKF8a+tW1VPPVVVRFMfeUTDwsL0zjvv9GjTnCu3vTEoXc7FYlu3bq3wvnL4ey6C008/Xc8555xi17/44osKVPhahkBz5MgRPemkk/SYY44psTptxYoV2rZtWw0ODtann3660MyRgcDfc7LksERh/Of771WbNnVGd509W19//XUFcruUlmbPnj1aq1YtrwwOd+KJJ+rZZ59d4f0U5M/ZzUaPHq1hYWFF1l+npKRoVFSU9u3b16918b7y888/a1hYmA4cOLDQ88vOztbXX39dw8PDtUWLFvrNN9/4KcrS+fvHRg5LFMY/3npLNSxM9fjjVdesUVXniuiOHTuW6Yvrqquu0qioqAo15q1Zs8br7R2BYObMmQroV199VWjdLbfcosHBwbrGfe2ro+eff14BnTx5cu6yAwcO5I5LdeGFF+rOnTv9GKFnAmEqVUsUpnJlZKiOGuV8fM4/X3X3blV1xr0B9IknnijT7pKSkhTQd999t9whPfzwwyoiusOL40UFgtTUVBURfeyxx/It//nnnzUoKEhvv/12P0VWObKysrR3794aGhqqzZs3VxHRkJAQBXTcuHEBWdUUqCqSKIIwpiz27nVmJXvhBbj9dmdWN3e6ypwpTePi4sq0y/POO4+YmBjefPPNcoWkqiQmJhIbG0uzZs3KtY9A1aBBA04++eR8U8qqKiNHjiQqKopHH33Uf8FVgqCgIAYNGkRGRgY7duxAVcnMzCQ8PJyYmBiCguwrrFKUN8P462ZnFH60bp1qhw6qISGqeaoCVJ0643bt2um5555brl0/9thjCpSru+CKFSsU0Ndff71cxw50t99+u0ZEROT2avroo48U0P/85z9+jqxyBEpjcFWHnVEYn1u0CM44w5lXeP58+Ne/8q1esmQJGzZs4J///Ge5dj98+HBEhLfffrvM2yYmJhIcHMxll11WrmMHutjYWA4fPsyyZcs4evQoo0ePpnPnztx8883+Dq1SbN26tUzLjfdZojCl27ABBgyAxo1h6VLo1atQkSlTplCrVi0uv/zych2idevWXHDBBbz11ltkZWV5vJ261U7nn38+jRs3LtexA10v9/VeuHAhEyZMYNOmTbz44ouEhIT4ObLK0bp16zItN95nicKUbNcu6N/fuf/ZZxATU6hIRkYG06dPZ+DAgdSvX7/ch7r++uvZunUrCxYs8HibZcuWsXnzZq688spyHzfQNWnShObNm/PII49w//33ExERwc6dO/0dVqUZP348kZGR+ZZFRkYyfvx4P0VU81iiMMU7fBgGDoTt22HWLGjbtshin332Gbt37+bqq6+u0OEuvfRSGjZsWKZG7cTEREJDQxk8eHCFjh3IEhIS2LlzJ+np6QAcPnyY+Ph4EhIS/BxZ5YiLi2Py5MlER0cjIkRHRzN58uQyd5owFVDexg1/3awxu5JkZaledpkzl/X775dY9IorrtDGjRsXO890WdyXBGhFAAAgAElEQVR+++0aFhamu90utyWHmKWtWrXSiy++uMLHDWTWmGu8AWvMNl53zz3wwQfw3HNQQiPx3r17mTlzJkOHDiU0NLTCh73hhhtIT0/n3XffLbXs4sWL2bZtW7WudgJrzDX+Z4nCFPbKK851ErfdBnfdVWLR999/n6NHj1a42inHKaecwmmnneZR9dP06dMJDw9n0KBBXjl2oLLGXONvHicKEekpIte595uISBvfhWX85pNP4M47nbaJCRNApMTiU6ZMoX379px++uleC+H6669n5cqVrFixotgyWVlZvPfee/Tv35969ep57diByBpzjb95lChE5BHgPuABd1EoMNVXQRk/WbIEhg2Drl1h2jQIDi6xeHJyMgsXLuTqq69GSkkoZXHVVVcRHh7O//73v2LLfPPNN/z555/VvtoJrDHX+J+nZxSDgYHAIQBV3QHU9VVQxg82b3aulTj2WKeHU4FfsEXJ6XVT3ovsitOgQQP+8Y9/kJCQwJEjR4osk5iYSGRkJJdccolXjx2o4uLi2LJlC9nZ2WzZssWShKlUniaKdLfVXAFEpLbvQjIFJSQk5I5rExMT4/1ukamp0K8fZGTAnDnQtGmpm6gqU6ZM4ZxzziGmiGsrKuqGG25gz549fPzxx4XWZWZm8sEHH3DJJZdQu7Z9FI3xNU8TxQwRmQREici/gCTgv74Ly+RISEggPj6e5ORkVJXk5GTv9qE/cgQuvdQ5o/jkEzjhBI82W758Ob/++qvXzyZynHfeeURHRxfZqP3ll1+SkpJSI6qdjAkEHiUKVX0OeB/4AOgAPKyq//FlYMYxZswY0tLS8i1LS0tjzJgxFd95djZcdx188w288w6cc47Hm06ZMoWwsDCuuOKKisdRhKCgIK677jrmz5/Pli1b8q2bPn06devWpV+/fj45tjEmP08bs9sA36jqPap6N/CtiMR4sN1FIrJORDaKyP3FlBkiIr+IyBoRKb3zfA3j0z70Y8bA9Onw1FMwdKjHm2VkZDBt2jQGDBhAgwYNKh5HMYYPHw7AO++8k7ssPT2dDz/8kEGDBhEREeGzYxtj/uZp1dN7QHaex1nusmKJSDAwEegHdAKGiUinAmXa4fSkOltVOwMjPYynxiiur7yqctNNN7F27dry7XjSJCdB3HQT3HtvmTb94osvSElJ8dq1E8WJjo6mT58+vPXWW2RnZ+cee+/evVbtZEwl8jRRhKhqes4D935YKducAWxU1U1u+elAwSuj/gVMVNU97n5rzkhnHiqqD32tWrU477zzeOedd+jUqRP9+/fniy++wOlv4IE5c2DECGewv1deKfVaiYKmTJlCo0aNKqXq5/rrryc5OTl3oMDExESioqLo27evz49tjHF4mihSRGRgzgMRGQTsKmWbFsC2PI+3u8vyag+0F5HvRGSxiFzkYTw1RlxcHHfeeSdAbh/6N954gwULFrBt2zYef/xxli9fTt++fTnllFN46623OHr0aPE7/PFHGDIETjkFEhOhjENV79+/n48//pgrr7ySsLDSfitU3KWXXkqDBg148803OXLkCB9//DGDBw+ulGMbYxyeJoqbgX+LyFYR2YZz8d1NXjh+CNAOOBcYBvxXRKIKFhKReBFZJiLLUlJSvHDYqic4OJh9+/bl60PfpEkTHnroIbZu3cpbb70FOL/Ao6OjGTt2LDmvVb7utaefTkJ4OMyeDXXqlDmODz74gCNHjvist1NBtWrVIi4ujo8++ohp06Zx4MABq3YyppJ52uvpN1U9E6etoaOqnqWqG0vZ7HegVZ7HLd1leW0HZqpqhqpuBtbjJI6Cx5+sqt1UtVuTJk08CblaSUpK4swzz6Ru3aKvcQwPD2f48OH89NNPfPHFF3Tt2pWHH36Y1q1b07t3b2688ca/u9dmZxOflkbCl1+WK5YpU6bQtm1bzjzzzIo8pTJp1qwZR48e5frrrycoKIi//vqr0o5tjMGzYcaBcOAq4N/Awzm3UrYJATYBbXDaM34COhcocxHwjnu/MU5VVaOS9lvThhlPTU1VEdFHHnmkTNv98ssvGh8fX+Tw1JRziOqtW7eqiOijjz5a5m3La+rUqRoZGZkv9sjISJ06dWqlxWBMdUAlDDP+CU5DdCbOMB45t5ISUCZwGzAPWAvMUNU1IvJ4nvaOecBuEfkF+BK4R1V3exhTjfDll1+iqvTp06dM23Xs2JFJkyYVOwZTcnIyN998My+//DJJSUns2LGj2MbwnKqr1q1bo6qVOgifT68jMcZ4RIr7cshXSGS1qp5YCfGUqlu3brps2TJ/h1FpRowYwZQpU0hNTS3XfA8xLVqQvGNHoeXh4eFERkayZ8+e3GX169enU6dOdOrUiY4dO9KpUyfWr1/Pv//973xf1pGRkZU2KF1QUFCRCUxEcrvMGmNKJyLLVbVbubb1MFFMBv6jqqvKcxBvqmmJon379rRv357Zs2eXfePMTBI6diR+40by/ibP+aK/6qqr2LlzJ7/88ku+29q1a0ttB4iOji50xbQvxMTEkJyc7LfjG1NdVCRReNpG8QuQDqwDfgZWAT+Xt76rIrea1EaxZcsWBfTFF18s3w4ee0wVdOqIERodHa0iotHR0R7V7+/atUu/+eabYts4RKR8MZWRtVEY4x1UoI3C00QRXdStvAetyK0mJYo333xTAV21alXZN160SDU4WDUurkIxBMJ8zVOnTi1zojPG5FeRROFp99hkVU0GDhf4wjA+lJSURNOmTencuXPZNjxwAOLioGVLmDixQjEEwuxqNheDMf7l6aCAA0VkA7AZ+BrYAsz1YVw1XnZ2NklJSfTp06fss8fdeSds2QJTp0L9+hWKw2ZXM8Z4On7DWOBMIElVTxWR84DKuTS3hlq9ejUpKSll7hbL++/DW285I8P27OmVWOLi4iwxGFODeXodRYY61zcEiUiQqn4JlK/13HgkKSkJgPPPP9/zjbZvh/h4OP10eOQRH0VmjKlpPD2j2CsidYCFQIKI7KSUC+5MxSQlJdGhQwdatWpVemFwJiG69lpIT4eEBCjHNRfGGFMUT88oBuE0ZN8FfAb8BgzwVVA1XXp6Ol9//XXZqp1eeAEWLIAJE6BdoeGyjDGm3Dw6o1DVvGcP7xRb0HjF4sWLSUtL8zxRrFwJ//43DB4MN9zg2+CMMTVOiYlCRL5V1Z4icoD83WEFUFWtvEF/apCkpCSCgoI499xzSy+clgZXXQWNG8N//1vmSYiMMaY0JSYKVe3p/i16fGvjE0lJSZx++ulERRWamqOwe++FtWvh88+hUSPfB2eMqXFKbaMQkWAR+bUygjGwb98+lixZ4lm105w5zgV1d90FF1zg++CMMTVSqYlCVbOAdSLSuhLiqfG+/vprsrKySk8UO3fCddfBSSfBE09UTnDGmBrJ0+6xDYA1IrKEPN1iVXVg8ZuY8khKSiIiIoIePXoUX0gVrr8e9u2D+fOhVq3KC9AYU+N4mige8mkUJldSUhK9evUiPDy8+EKvvw6ffgovvQQnBsQ0IcaYaszT7rFf+zoQA7///jtr167l+uuvL77Q2rUwejRceCHcfnvlBWeMqbE8HRTwTBFZKiIHRSRdRLJEZL+vg6tp5s+fD1B8+0R6ujMqbO3aznhO1hXWGFMJPK16egUYCryHM8bTNUB7XwVVUyUlJdG4cWNOPvnkogs89BCsWAGffALNmlVucMaYGsvTITxQ1Y1AsKpmqepbwEW+C6vmUVWSkpI4//zzCQoq4m1ZuBCefdYZ9G+g9SEwxlQeT88o0kQkDFgpIs8Af1CGJGNKt3btWv7444+iq52OHnUSREyMM6aTMcZUIk+/7K92y96G0z22FXCZr4KqiXKGFb+gqAvnnn0W1q2DV1912ieMMaYSeXpG0RX4VFX3A4/5MJ4aKykpibZt2xIdHZ1/xcaNMG4cDBkCF1ltnzGm8nl6RjEAWC8iU0TkEhHxNMEYD2RkZPDVV18VrnZShVtvhbAwePFF/wRnjKnxPEoUqnod0Ban19Mw4DcRecOXgdUkS5cu5cCBA4UTxXvvOYP9jR8PzZv7JzhjTI3n8ZmBqmaIyFyc4cYjgEuBG30VWE2SlJSEiHDeeef9vXDfPrjzTujaFUaM8F9wxpgaz9ML7vqJyNvABpxG7DeAY30YV42SlJRE165dadiw4d8LH3zQGfhv0iQIDvZfcMaYGs/TNoprgI+BDqo6XFXnqGqmD+OqMQ4ePMiiRYvyVzstW+YMH37rrc4ZhTHG+JGnYz0N83UgNdXChQvJzMz8O1FkZcFNN8Gxx8LYsf4NzhhjKPtUqJL3r02FWnFJSUnUqlWLs88+21kwcSL8+CMkJkL9+v4NzhhjsKlQ/S4pKYmePXtSq1Yt+P13p23iwgvhiiv8HZoxxgBl6PUkIg1wrsjO3UZVf/RFUDXFn3/+yapVq4iLi3MW3HUXZGQ4ZxU2MqwxJkB4lChEZCwwHNgEZLuLFejtm7BqhgULFgDusOJz5zrXTYwbB8cf7+fIjDHmb56eUQwBjlfV9LLsXEQuAl4CgoE3VPWpYspdBrwPnK6qy8pyjKosKSmJhg0b0qVDB6eq6YQT4O67/R2WMcbk42miWA1EATs93bGIBAMTgQuA7cBSEZmpqr8UKFcXuBP4wdN9Vwc5w4r37t2b4CefhM2b4csvoaQpUI0xxg88vY7iSWCFiMwTkZk5t1K2OQPYqKqb3DOR6cCgIsqNBZ4GjngcdTWwYcMGtm3bRp/OnZ3RYa+9Fs49199hGWNMIZ6eUbyD82W+ir/bKErTAtiW5/F2oHveAiJyGtBKVT8VkXuK25GIxAPxAK1bt/bw8IEtZ1jxPnPmQJ06TrIwxpgA5PHERar6sjcPLCJBwAs4jeQlUtXJwGSAbt26qTfj8JekpCRiGjfmuKVL4b//hSZN/B2SMcYUydNE8Y2IPAnMBI7mLCyle+zvON1pc7R0l+WoC5wIfCVOV9BjgZkiMrC6N2hnZWWxYP58rjh6FDnrLLj+en+HZIwxxfI0UZzq/j0zz7LSuscuBdqJSBucBDEUuCp3Y9V9QOOcxyLyFXB3dU8SAMuXL2ff/v30CQqC11+HoubINsaYAOHpWE/nlV6q0DaZInIbMA+ne+z/VHWNiDwOLFPV0hrDq62kN5ypPHqPGAEnneTnaIwxpmSiWnqVv4g0BZ4AmqtqPxHpBPRQ1Td9HWBB3bp102XLqvBJR0YGvRs0YE96Oiv27LE5sI0xlUJElqtqt/Js62mdx9s4ZwY506ytB0aW54A1XdpTT/HdoUP0ueQSSxLGmCrB00TRWFVn4HaNdeeiyPJZVNXVli18O24c6UCfm27ydzTGGOMRTxPFIRFphNOAjYicCezzWVTVkSrcdhtJqoSFhdGzZ09/R2SMMR7xtNfTKJyusceLyHdAE+Byn0VVHX3yCXz6KUktWnBWu3bUtmonY0wVUeIZhYicLiLHutdLxAL/xrmO4nOcK62NJ44ehdGj2dWhAyt+/z3/tKfGGBPgSqt6mgTkjBh7FjAGZ6C/PbhXShsP/Oc/sGkTCy67DMAShTGmSimt6ilYVVPd+1cCk1X1A+ADEVnp29CqiZQUZ+7r/v1JSkmhfv36dO3a1d9RGWOMx0o7owgWkZxkcj6wIM86j2fHq9EeeYSEAweIWbGC//73v6Snp5OYmOjvqIwxxmOlfdlPA74WkV3AYeAbABFpi/V6Kt2aNSS89hrxISGk/fEHAIcPHyY+Ph7g7ylQjTEmgJV6ZbbbFbYZ8LmqHnKXtQfq+GPO7Cp1ZfZFFxHz+eckF/EaR0dHs2XLlsqPyRhTI1XkyuxSq49UdXERy9aX52A1yty5MG8eW4tZvXVrcWuMMSaw2LClvpCRAaNHo23bUqdu3SKLVJcJmIwx1Z8lCl+YPBldu5aRJ5zAgQMHCAnJf+IWGRnJ+PHj/RScMcaUjSUKb9uzB334YUa2bMnLs2czatQo3n77baKjoxERoqOjmTx5sjVkG2OqDOvi6mU6dix3pqbyn9RURo8ezbPPPouIWGIwxlRZdkbhRbp+PXe89BL/Ae6+++7cJGGMMVWZnVF4iapy+0UXMTE7m7tvuYVnnnnGkoQxplqwMwovUFVuHzyYiZs3c8855/DMxImWJIwx1YYligpSVW679VYmfvIJ99aty9OffWZJwhhTrViiqIDs7GxuvfVWXn3tNe4Fnpo0CYmM9HdYxhjjVdZGUU7Z2dncdtttvPbaa9xXuzZPnnQSMnSov8MyxhivszOKcsg5k3jttde4v0cPnjx0CJkwAazKyRhTDVmiKKPs7GxGjBjB66+/zv233MITy5cjV10F3bv7OzRjjPEJSxQeSEhIICYmhqCgIOrXr8+kSZN44IEHeCI1FQkKgief9HeIxhjjM9ZGUYqEhATi4+NJS0sD4ODBg4SEhNA5LAxJTIQHHwQb4M8YU42VOh9FoKns+ShiYmJITk4utDw6LIwtjRrB+vVQp06lxWOMMeVRkfkorOqpFMXNG7E1PR2eeMKShDGm2rNEUYri5o1oHRoK11xTydEYY0zls0RRivHjxxMaGppvWSQw/r77IMhePmNM9WffdKWIi4ujZcuWhIWFOfNJiDC5Wzfixo71d2jGGFMpLFGUYtOmTWzevJlx48aRPXw4W0JDiZs+3d9hGWNMpbFEUYr33nsPgCGdOsHbb8Mdd8Dxx/s3KGOMqUQ+TRQicpGIrBORjSJyfxHrR4nILyLys4jMF5FoX8ZTHomJiXTv3p3oCROgYUMYM8bfIRljTKXyWaIQkWBgItAP6AQME5FOBYqtALqp6snA+8AzvoqnPDZs2MCKFSu4sksXSEpyLq6LivJ3WMYYU6l8eUZxBrBRVTepajowHRiUt4Cqfqmqae7DxUBLH8ZTZomJiQBc8e23EBMDt9zi34CMMcYPfDmERwtgW57H24GSRs67AZhb1AoRiQfiofjrGnxhxowZ9OzQgZZr1sDUqRAeXmnHNsaYQBEQjdki8k+gG/BsUetVdbKqdlPVbk2aNKmUmNauXcuqVasYkpICp5wCw4ZVynGNMSbQ+PKM4negVZ7HLd1l+YhIH2AMEKuqR30YT5kkJiYiIlyemgrvvmsX1xljaixffvstBdqJSBsRCQOGAjPzFhCRU4FJwEBV3enDWMpEVUmcNo3YkBCa9e4Nffv6OyRjjPEbnyUKVc0EbgPmAWuBGaq6RkQeF5GBbrFngTrAeyKyUkRmFrO7SrV69Wp+Xb+eKzMy4KmnbOY6Y0yN5tP5KFR1DjCnwLKH89zv48vjl1fim28SBPxj4EA4/XR/h2OMMX5lFe8FqCqJb79NbxGOee45f4djjDF+Z4migBUff8zGffu4MjYW2rXzdzjGGON3ligKmPHAA4QAgydO9HcoxhgTECxR5KE//EDiunX0Of54GnUqONqIMcbUTJYocqiydMQItgBX3n23v6MxxpiAYYkix7x5JP74I2EhIVw6dKi/ozHGmIBhiQIgO5vse+/lveBgLrzwQqJshFhjjMlliQIgIYHFq1axLSuLIXY2YYwx+ViiOHIEHnyQxCZNCA8PZ+DAgaVvY4wxNYglitdeI2vrVt7LyqJ///7Uq1fP3xEZY0xAqdmJYt8+GDeO77p144/UVK688kp/R2SMMQGnZieKZ56B1FQS27QhIiKCiy++2N8RGWNMwKm5iWLHDnjxRTKvvJL3v/6aSy65hDp16vg7KmOMCTg1N1E89hhkZvL1xRezc+dOq3Yyxphi1MxE8euv8OabcMstzPj2W2rXrk3//v39HZUxxgQkn85HEbDGjIHISDLuu48PTj6ZgQMHEhER4e+ojDEmINW8M4pFi+DDD+Gee1iwahW7d++2aidjjClBzUoUqnDffdC0Kdx1F4mJidSrV4+LLrrI35EZY0zAqllVT59+Ct98A6++SnpYGB999BGXXnop4eHh/o7MGGMCVs05o8jKgvvvh7Zt4cYb+eKLL9i7dy9Dhgzxd2TGGBPQas4ZxZQpsGYNzJgBoaEkJibSoEEDLrjgAn9HZowxAa3mnFG0aQPXXQeXX86RI0f4+OOPGTx4MGFhYf6OzBhjAlrNOaOIjXVuwLx58zhw4ID1djLGGA/UnDOKPBITE2nUqBHnnXeev0MxxpiAV+MSRVpaGjNnzuSyyy4jNDTU3+EYY0zAq3GJYu7cuRw6dMiqnYwxxkM1LlEkJiZyzDHHEOu2VxhjjClZjUoUBw8eZPbs2Vx++eUEBwf7OxxjjKkSalSimD17NocPH7ZqJ2OMKYMakSgSEhKIiYlh2LBhBAcHs3XrVn+HZIwxVUa1v44iISGB+Ph40tLSAMjKyuKmm25CRIiLi/NzdMYYE/iq/RnFmDFjcpNEjrS0NMaMGeOniIwxpmqp9omiuGomq34yxhjP+DRRiMhFIrJORDaKyP1FrA8XkUR3/Q8iEuPtGFq3bl2m5cYYY/LzWaIQkWBgItAP6AQME5FOBYrdAOxR1bbAi8DT3o5j/PjxREZG5lsWGRnJ+PHjvX0oY4yplnx5RnEGsFFVN6lqOjAdGFSgzCDgHff++8D5IiLeDCIuLo7JkycTHR2NiBAdHc3kyZOtIdsYYzzky15PLYBteR5vB7oXV0ZVM0VkH9AI2JW3kIjEA/FQviqjuLg4SwzGGFNOVaIxW1Unq2o3Ve3WpEkTf4djjDE1ii8Txe9AqzyPW7rLiiwjIiFAfWC3D2MyxhhTRr5MFEuBdiLSRkTCgKHAzAJlZgLXuvcvBxaoqvowJmOMMWXkszYKt83hNmAeEAz8T1XXiMjjwDJVnQm8CUwRkY1AKk4yMcYYE0B8OoSHqs4B5hRY9nCe+0eAK3wZgzHGmIqRqlbTIyIpQHIxqxtToMdUDWPPv2Y/f7DXwJ5/8c8/WlXL1RuoyiWKkojIMlXt5u84/MWef81+/mCvgT1/3zz/KtE91hhjjP9YojDGGFOi6pYoJvs7AD+z529q+mtgz98HqlUbhTHGGO+rbmcUxhhjvMwShTHGmBJVm0RR2iRJVZWItBKRL0XkFxFZIyJ3ussbisgXIrLB/dvAXS4i8rL7OvwsIqfl2de1bvkNInJtcccMNCISLCIrRGS2+7iNO9HVRnfiqzB3ebETYYnIA+7ydSJyoX+eSfmISJSIvC8iv4rIWhHpUcPe/7vcz/5qEZkmIrWq82dARP4nIjtFZHWeZV57v0Wkq4iscrd5WcSDqR1UtcrfcIYI+Q04DggDfgI6+TsuLz23ZsBp7v26wHqciaCeAe53l98PPO3e7w/MBQQ4E/jBXd4Q2OT+beDeb+Dv5+fhazAKeBeY7T6eAQx1778O3OLeHwG87t4fCiS69zu5n4lwoI37WQn29/Mqw/N/B7jRvR8GRNWU9x9nKoLNQESe9354df4MAL2A04DVeZZ57f0Glrhlxd22X6kx+ftF8dIL2wOYl+fxA8AD/o7LR8/1E+ACYB3QzF3WDFjn3p8EDMtTfp27fhgwKc/yfOUC9YYz6vB8oDcw2/1w7wJCCr73OOOK9XDvh7jlpODnIW+5QL/hjKi8GbfjScH3tQa8/zlz1jR039PZwIXV/TMAxBRIFF55v911v+ZZnq9ccbfqUvVU1CRJLfwUi8+4p9GnAj8ATVX1D3fVn0BT935xr0VVfY0mAPcC2e7jRsBeVc10H+d9HvkmwgJyJsKqqs8dnF+/KcBbbvXbGyJSmxry/qvq78BzwFbgD5z3dDk16zMA3nu/W7j3Cy4vUXVJFNWeiNQBPgBGqur+vOvU+WlQ7fo5i8glwE5VXe7vWPwoBKca4jVVPRU4hFP1kKu6vv8Abl38IJyE2RyoDVzk16D8zB/vd3VJFJ5MklRliUgoTpJIUNUP3cV/iUgzd30zYKe7vLjXoiq+RmcDA0VkC86c672Bl4AocSa6gvzPo7iJsKric8+xHdiuqj+4j9/HSRw14f0H6ANsVtUUVc0APsT5XNSkzwB47/3+3b1fcHmJqkui8GSSpCrJ7ZHwJrBWVV/IsyrvpE/X4rRd5Cy/xu0NcSawzz1lnQf0FZEG7q+0vu6ygKWqD6hqS1WNwXlPF6hqHPAlzkRXUPi5FzUR1kxgqNsjpg3QDqdBL+Cp6p/ANhHp4C46H/iFGvD+u7YCZ4pIpPu/kPP8a8xnwOWV99tdt19EznRfz2vy7Kt4/m608WLjT3+cHkG/AWP8HY8Xn1dPnNPMn4GV7q0/Tr3rfGADkAQ0dMsLMNF9HVYB3fLs63pgo3u7zt/PrYyvw7n83evpOJx/8o3Ae0C4u7yW+3iju/64PNuPcV+TdXjQyyOQbkAXYJn7GfgYpxdLjXn/gceAX4HVwBScnkvV9jMATMNpj8nAOaO8wZvvN9DNfS1/A16hQEeJom42hIcxxpgSVZeqJ2OMMT5iicIYY0yJLFEYY4wpkSUKY4wxJbJEYYwxpkSWKEy1IiIHfbDPLe5om6vEGcV3nIjUctc1F5H3S9g2SkRGeDsmYyqTdY811YqIHFTVOl7e5xac/um73KFUJgMZqlrqUN3u+FyzVfVEb8ZkTGWyMwpT7YnIAHdughUikiQiTd3lTdyx/de4g+0li0jjkvalqgeBm4FL3TkCYnLmDRCRziKyRERWunMDtAOeAo53lz0rInVEZL6I/OieoQxyt40RZ66J/7rxfC4iEe66tm7cP7nbHe8uv0dElrrHesx3r6Cp6SxRmJrgW+BMdQbVm44zGi3AIzhDPHTGGUOptSc7U2dQxs04w0DkdTPwkqp2wbn6dTvOAH6/qWoXVb0HOAIMVtXTgPOA5/NMHNMOmOjGsxe4zF2e4C4/BTgL+ENE+rrlz8C5cruriPTy+BUxpgxCSi9iTJXXEkh0B1MLw/mSB2d4lMEAqvqZiOwpwz6LmhVsETBGRFoCH6rqhiImDxPgCfdLPRtniOecIaM3q+pK9zZEEi8AAAGUSURBVP5yIEZE6gItVPUjN84jAG6i6AuscMvXwUkcC8vwHIzxiJ1RmJrgP8ArqnoScBPOeEDl5n55x+CMLZZLVd8FBgKHgTki0ruIzeOAJkBX98zjrzzxHM1TLouSf8gJ8KR7ptJFVduq6pvleT7GlMYShakJ6vP3UMp5G6C/A4ZA7i/0BqXtyG3MfhX4WFX3FFh3HLBJVV/GGZHzZOAAzhS2eWPZqaoZInIeEF3S8VT1ALBdRC51jxEuIpE4o4Ne78aDiLQQkWNKi9+Y8rCqJ1PdRIpI3hm8XgAeBd5zq5YW4EyCA86opNNE5GqcaqM/cb7Yi/Kl25YQBHwEjC2izBDgahHJcPf1hKqmish3boP3XOBpYJaIrMIZEfZXD57T1cAkEXkcZ0TRK1T1cxHpCCxyq7cOAv/k73kKjPEa6x5raiwRCQeyVDVTRHrgzCLXxd9xGRNo7IzC1GStgRkiEgSkA//yczzGBCQ7ozDGGFMia8w2xhhTIksUxhhjSmSJwhhjTIksURhjjCmRJQpjjDEl+n+OM7CzoPlm5AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "svm = model.semivariance(model.spherical, (4000, sill))\n",
    "geoplot.semivariogram(P, lags, tolerance, model=svm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The covariance modeling function function will return a spherical covariance model that takes a distance as input, and returns an covariance estimate. We've used the global variance of the porosity in `ZoneA.dat` as the sill."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "covfct = model.covariance(model.spherical, (4000, sill))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then krige the data, using the covariance model, the point we are interested in, (2000,47000), and `N=6` signifying that we only want to use the six nearest points. The output of the simple and ordinary kriging functions below is the krigin estimate, and the standard deviation of the kriging estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(12.826498066190497, 0.49628540680221894)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kriging.simple(P, covfct, pt, N=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(12.931216562350736, 0.4986494585776898)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kriging.ordinary(P, covfct, pt, N=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "est, kstd = kriging.krige(P, covfct, [[2000,4700],[2100,4700],[2000,4800],[2100,4800]], 'simple', N=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[12.82649807],\n",
       "       [12.71773246],\n",
       "       [12.83938422],\n",
       "       [12.73159357]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.49628541],\n",
       "       [0.48591652],\n",
       "       [0.49470641],\n",
       "       [0.48637714]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kstd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
